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ABSTRACT 


A  time  dependent,  three  dimensional  finite  element  approach  to  the  develop¬ 
ment  of  a  perfectly  matched  layer  for  numerical  calculations  of  surface  wave  radi¬ 
ation  in  a  half  space  is  presented.  The  development  of  this  new  element  required 
the  coupling  of  a  system  of  linear,  second-order,  partial  differential  equations  which 
describe  elastic  wave  propagation,  together  with  their  related  boundary  conditions, 
into  a  single  weak-form  (Galerkin)  wave  equation,  from  which  the  characteristics  of 
a  composite  finite  element  matching  layer  were  derived.  An  important  problem  of 
interest,  and  the  motivation  for  this  work,  is  the  optimization  of  a  source  for  use  in 
a  seismo-acoustic  sonar  for  the  detection  of  buried  mines.  Validation  of  the  perfectly 
matched  layer  occurs  by  employing  it  in  a  finite  element  analysis  to  compute  the  ra¬ 
diation  from  a  particular  transient  seismo-acoustic  source  array  and  showing  that  the 
results  agreed  with  the  results  of  previous  field  experiments  using  the  same  source 
performed  by  Naval  Postgraduate  School  students.  Various  source  excitations  are 
presented  which  maximize  the  energy  of  the  unidirectional  Rayleigh  wave  while  sup¬ 
pressing  the  energy  of  associated  body  waves.  Radiation  characteristics  are  analyzed 
in  a  linear,  isotropic,  homogeneous  half  space  with  a  discrete  number  of  transient 
seismic  sources.  The  hp-adaptive  finite  element  code  SAFE-T  (Solid  Adaptive  Finite 
Element  -  Transient),  a  Finite  Element  Method  (FEM)  implementation  developed 
by  the  author  utilizing  Altair  Engineering’s  Prophlcx  kernel,  is  used  to  perform  the 
numerical  computations.  Results  for  radial  and  vertical  wave  strengths  are  given  in 
terms  of  their  total  displacement  magnitudes.  This  work  represents  an  important  step 
forward  in  the  development  of  tools  needed  to  pursue  seismo-acoustic  sonar  technol¬ 
ogy  for  buried  mine  detection,  as  well  as  for  the  analysis  of  all  three-dimensional, 
time-dependent  elasto- dynamic  problems. 
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I. 


INTRODUCTION 


A.  BACKGROUND 

The  original  analytic  development  of  a  surface  displacement  resulting  from  a 
vertical  surface  disturbance  is  clue  to  Lamb  [Ref.  1],  His  method  intricately  syn¬ 
thesized  the  solution  of  a  point  and  line  pressure  pulse  varying  with  nearly  impulse 
time  dependence  from  the  periodic  solution,  and  yielded  displacements  of  the  surface 
excited  by  a  transient  unit  step  source.  Lamb’s  work  has  been  extended  and  explored 
by  many  authors;  Pekeris(1955)  [Ref.  2],  Garvin(1960)  [Ref.  3],  and  Graff(1975)  [Ref. 
4]  just  to  name  a  few.  J.  D.  Achenbach  [Ref.  5]  notes  in  his  well  known  book  Wave 
Propagation  in  Elastic  Solids,  “In  recent  years  the  methods  and  solutions  in  Lamb’s 
paper  have  been  cast  in  a  somewhat  more  elegant  form  and  more  detailed  compu¬ 
tations  have  been  carried  out,  particularly  for  loads  of  arbitrary  time  dependence”. 
Of  particular  note  are  the  solutions  worked  out  by  Pekeris  [Ref.  2],  He  explored 
closed-form  solutions  for  vertical  and  horizontal  surface  displacements  in  response  to 
a  transient  vertical  point  surface  pressure.  However,  the  solutions  were  only  for  unit 
step  and  delta  functions  which  are  not  physically  realistic,  but  rather  represented 
limiting  cases.  H.  M.  Mooney  [Ref.  6]  demonstrated  the  effects  of  changing  Poisson’s 
ratio  in  the  domain,  and  showed  how  to  analyze  transient  arbitrary  sources.  Mooney 
presented  the  results  in  a  general  form  which  permit  convenient  application  to  other 
problems,  and  most  importantly  comparisons  to  numerical  methods  as  well.  The 


Figure  1.  Surface  Wave. 
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Finite  Element  Method  (FEM)  is  a  popular  tool  for  the  numerical  approximation 
of  partial  differential  equations  (PDEs).  Because  of  the  availability  of  powerful  and 
inexpensive  computing  platforms,  the  held  of  computational  mechanics  has  come  to 
rely  greatly  on  FE  methods  to  numerically  solve  PDEs  which  arise  in  the  study  of 
various  disciplines  within  the  physical  sciences  [Ref.  7].  Of  particular  interest  to  the 
problem  at  hand  is  wave  motion  in  an  elastic  medium.  Seismic  wave  motions  occur  in 
solid  media  due  to  clastic  forces  present  in  and  around  solids.  A  main  peculiarity  of 
elastic  seismic  waves  is  that  in  an  isotropic  medium  there  can  be  two  or  more  types  of 
waves  that  may  travel  with  different  velocities  and  different  polarizations.  This  is  clue 
to  different  clastic  states  of  deformation:  mainly  shear  and  compression.  Seismic  or 
clastic  techniques  have  shown  considerable  promise  in  the  reliable  detection  of  various 
types  of  buried  objects  [Ref.  8].  In  particular,  the  study  of  seismic  source  array  sys¬ 
tems  to  generate  surface  waves  for  the  detection  of  land  mines  in  an  clastic  half-space 
has  provided  a  unique  application  for  matching  mathematical  FEM  models  to  results 
derived  through  direct  experimentation  [Ref.  9].  The  end  goal  of  the  mathematical 
model  is  two-fold.  First,  we  want  to  determine  an  optimal  position,  i.e.,  spacing,  of 
seismic  transient  sources,  and  second,  we  want  to  accurately  time  the  excitation  of 
seismic  sources  as  to  profit  from  the  destructive  and  constructive  interference  that 
results. 


B.  COMPUTER  MODELS 

A  primary  obstacle  in  formulating  an  accurate  FEM  half-space  model  is  ad¬ 
dressing  how  to  simulate  wave  phenomena  for  the  entire  unbounded  clastic  medium. 
Computer  simulations  are  finite  so  significant  steps  must  be  taken  to  truncate  these 
systems  that  attempt  to  model  infinite  or  semi-infinite  domains.  Computers  are  finite 
machines  with  limited  resources  i.e.,  memory,  hard  drive  space,  cpu  speed  etc.  Mod¬ 
eling  an  infinite  domain  on  a  computer  system  is  like  attempting  to  model  the  effects 


2 


Contour  (Analysis  system) 
Displacement  DO 
■— 8.269E-02 
1-7.000E-02 
—  6UU0E-02 
I — 5.000E-02 
S-4.000E-02 
3.000E-02 

- 2.000E-02 

M— 1.000E-02 
■— 0.000E*00 
®--1.073E*01 


Figure  2.  Axisymmetric  wedge  of  a  three  dimensional  half-space.  Modeling 
an  infinite  domain  on  a  computer  system  is  like  attempting  to  model  the 
effects  and  properties  of  the  ocean  in  a  bucket.  Axisymmetry  is  a  technique 
used  to  conserve  resources  while  not  diminishing  the  scope  of  the  problem. 

and  properties  of  the  ocean  in  a  bucket.  The  challenge,  then,  is  to  use  a  finite  space 
to  accurately  model  an  infinite  one.  Several  approaches  have  been  used  to  accomplish 
this  task.  Figure  2  is  an  attempt  to  model  wave  propagation  by  taking  advantage 
of  the  axisymmetric  properties  of  the  problem.  Yet,  without  significant  computer 
resources  to  model  an  enormous  domain,  even  this  method  falls  short  of  accurately 
modeling  an  infinite  half-space.  In  their  paper,  “Absorbing  Boundary  Conditions  For 
Acoustic  And  Elastic  Wave  Equations,”  Clayton  and  Engquist(1977)  [Ref.  12]  de¬ 
veloped  an  absorbing  boundary  condition  for  the  2D  elastic  wave  equation  using  a 
paraxial  approximation  method.  With  an  absorbing  boundary,  the  domain  need  not 
be  as  large  and  a  solution  moves  closer  to  within  reach.  Research  continues,  how¬ 
ever,  in  improving  and  extending  numerical  methodologies  of  terminating  elastic  and 
acoustic  unbounded  domains.  The  challenge  has  led  many  to  attempt  the  use  of 
traditional  frequency  domain  models  in  efforts  to  obtain  time  domain  solutions.  One 
clever  approach  was  explored  by  Bernal  and  Youssef(1998)  [Ref.  13]  who  improved 
the  use  of  a  hybrid  frequency-time  domain  procedure  which  iterates  between  the  fre¬ 
quency  and  time  domain.  A  reference  linear  system  is  solved  in  the  frequency  domain 
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while  the  equations  of  motion  are  solved  in  the  time  domain  with  the  unbounded 
medium  represented  by  frequency  independent  springs  and  dampers.  The  frequency 
dependency  of  the  impedance  of  these  springs  and  dampers  is  introduced  into  the 
system  by  means  of  frequency  domain  force  evaluations  at  the  end  of  each  time  step. 
Basu(2003)  [Ref.  14]  comments  that  this  method  is  computationally  demanding,  and 
requires  careful  implementation  to  ensure  stability. 


C.  THE  PERFECTLY  MATCHED  LAYER 

A  modern  silent  boundary  condition  that  is  receiving  a  flurry  of  attention  is 
the  perfectly  matched  layer  method.  First  introduced  by  Berenger  [Ref.  15]  in  1994  for 
use  with  electromagnetic  waves,  and  almost  immediately  applied  by  Chew  and  Wee- 
ton[Ref.  16]  for  use  with  Maxwell’s  equations,  it  has  been  shown  to  absorb  completely 
incident  plane  waves  without  reflection  over  a  broad  range  of  incidence  angles  and 
frequencies.  This  method  is  growing  rapidly  and  has  been  used  in  many  fields,  rang¬ 
ing  from  use  with  eddy-current  problems  by  Kosmanis’(1999)  [Ref.  17]  to  application 
to  electromagnetic  media  by  Cummer’s(2003)  [Ref.  18]  and  linearized  shallow  water 
equation  models  by  Navon,  Neta,  and  Hussaini  (2001)  [Ref.  19].  The  reason  is  that 
the  PML  can  be  formulated  at  a  relatively  small  computational  cost.  Recently,  Basil 
and  Chopra (2003)  [Ref.  11]  developed  the  PML  concept  in  terms  of  time-harmonic 
elastodynamics  by  utilizing  insights  from  electromagnetics  and  presented  a  novel  dis¬ 
placement  based  finite  element  (FE)  implementation  for  time-harmonic  plane  strain 
and  three-dimensional  analysis.  Later  Basil  and  Chopra(2004)  [Ref.  14]  extended  the 
idea  to  a  2D  transient,  displacement-based,  finite  element  (FE)  method  implemen¬ 
tation  for  anti-plane  and  for  plane-strain  motion  by  a  special  choice  of  a  coordinate 
stretching  technique. 

The  basic  idea  behind  the  PML  methodology  is  to  surround  the  computational 
domain  at  the  infinite  boundary  with  some  type  of  absorbing  layer.  As  you  can  see 
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Figure  3.  Perfectly  matched  layers  of  a  computational  domain.  The  com¬ 
putational  domain  is  in  the  center  surrounded  by  a  boot-like  anisotropic 
sponge  which  is  perfectly  matched  at  the  interface  between  the  two  do¬ 
mains. 


from  Figure  3,  the  boundary  layer  is  composed  of  the  same  kind  of  elements  as  the 
original  computational  domain.  The  formulation  of  the  elements  is  the  same  for  both 
the  computational  domain  and  the  absorbing  boundary,  but  the  boundary  layer  has 
slightly  different  properties.  This  boundary  layer  will  be  referred  to  as  a  perfectly 
matched  layer  (PML)  which  is  in  substance  a  perfectly  matched  media  (PMM). 


D.  DISSERTATION  GOAL 

It  is  well  established  that  analytic  procedures  [Ref.  10]  and  their  software 
implementations  incorporating  surface  wave  interactions  are  currently  formulated  in 
the  frequency  domain  for  three  dimensions  and  for  up  to  two  dimensions  in  the  time 
domain  [Ref.  11],  This  dissertation  presents  a  three  dimensional  time  dependent 
finite  element  approach  to  the  development  of  a  perfectly  matched  layer  for  numerical 
calculations  of  surface  wave  radiation  in  a  half  space.  Various  source  excitations 
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will  be  examined  to  maximize  the  energy  of  the  unidirectional  Rayleigh  wave  while 
suppressing  the  energy  of  associated  body  waves.  The  radiation  characteristics  are 
analyzed  in  a  linear,  isotropic,  homogeneous  half-space  with  a  discrete  number  of 
transient  seismic  sources.  The  mathematical  formulation  of  the  problem  consists  of 
the  coupling  of  a  system  of  linear,  second  order,  partial  differential  equations  and 
related  boundary  conditions  into  one  single  wave  equation  from  which  a  composite 
clastic  element  is  derived.  A  time  dependent  3D  perfectly  matched  layer  (PML) 
is  developed  to  handle  the  semi-infinite  half  space.  Outgoing  waves  are  completely 
absorbed  in  the  PML  with  minimal  reflections  from  all  angles  of  incidence.  The  Finite 
Element  Method  (FEM)  using  the  hp-adaptive  finite  element  kernel,  ProPHLEX,  is 
used  to  perform  the  numerical  computations. 

The  following  is  a  brief  summary  of  the  remainder  of  this  dissertation.  In 
Chapter  II,  the  major  steps  involved  with  developing  both  analytic  and  numerical 
solutions  to  clastic  wave  motion  problems  are  outlined.  Beginning  with  the  derivation 
of  the  elastic  partial  differential  equations  used  to  model  wave  motion,  methods  of 
solving  the  elastic  PDE  will  focus  primarily  on  solutions  which  allow  the  forcing 
function  or  source  to  be  arbitrary.  A  brief  introduction  to  the  finite  element  method 
will  be  presented  where  boundary  and  initial  conditions  are  constrained  to  produce 
a  well-posed  and  useful  mathematical  model.  Chapter  II  ends  with  a  discussion  of 
the  time  marching  scheme  and  the  mesh  conditions  employed  to  make  the  model 
dynamic.  Chapter  III  focuses  primarily  on  the  complex  problem  of  truncating  the 
computational  domain  to  only  allow  minimal  reflections  from  the  boundary.  A  terse 
review  of  the  methodology  involved  in  the  implementation  of  PMLs  in  the  frequency 
domain  is  given  followed  by  a  more  robust  three  dimensional  transformation  of  the 
vector  quantities  to  the  time  domain.  Considerable  time  is  spent  building  the  weak 
or  Galerkin  forms  of  the  equations.  This  is  necessary  in  order  to  implement  the 
mathematical  model  into  the  finite  element  code.  Chapter  III  concludes  with  SAFE- 
T’s  stunning  results  which  showcase  the  effectiveness  of  the  transient  PML  and  some 
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sensitivity  analysis.  Chapter  IV  is  focussed  on  the  specific  application  of  the  end-fire 
array.  SAFE-T  is  used  to  analyze  characteristics  of  an  end-fire  array  to  determine  the 
optimum  space  and  timing  that  results  in  the  greatest  magnitude  of  the  Rayleigh  wave 
while  suppressing  other  surface  and  non-surface  waves.  Finding  and  conclusions  are 
presented  in  Chapter  V  along  with  a  discussion  about  possible  future  contributions 
based  on  the  findings  of  this  investigation. 
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II. 


WAVE  MOTIONS  IN  ELASTIC  MEDIA 


A.  INTRODUCTION 

This  chapter  presents  the  results  of  a  newly  developed  time  dependent  finite 
element  (FE)  computer  code  capable  of  solving  3D  vector-valued  clastic  partial  differ¬ 
ential  equations.  A  review  of  techniques  used  to  find  analytic  solutions  is  discussed  for 
the  surface  response  to  an  instantaneous  transient  point  load  located  one  meter  away 
from  the  source.  Finding  analytic  closed-form  solutions  for  arbitrary  surface  traction 
are  rare.  The  analysis  of  the  response  to  a  point  or  line  load  acting  on  a  previously 
tranquil  half-space  is  more  common.  Care  is  taken  in  this  chapter  to  include  some 
of  the  difficulties  involved  in  forming  point  or  line  loaded  elastic  partial  differential 
equation  models  while  at  the  same  time  highlighting  the  various  mathematical  tools 
available  to  overcome  those  difficulties.  Integral  transformations,  most  useful  when 
an  appropriate  inversion  can  be  found  or  evaluated,  will  play  a  role  in  simplifying 
the  model  and  hireling  the  correct  solution.  One  major  challenge  to  the  point  load 
problem,  however,  is  that  the  solution  is  non-physical.  The  point  source  is,  in  fact,  a 
singularity  which  propagates  along  the  free  surface  of  the  mathematical  model.  Sin¬ 
gularities  are  extremely  difficult,  if  not  impossible,  to  model  numerically  because  they 
involve  quantities  which  are  not  finite.  An  integral  convolution  is  utilized  to  achieve 
a  non-singular  result  by  discretely  convolving  an  arbitrary  transient  source  with  the 
non-physical  solution  to  the  point  load  problem.  The  result  is  a  new  analytic  solution 
that  is  both  realistic  and  useful.  The  displacements  from  the  discrete  convolution  will 
be  compared  to  the  response  calculated  by  the  FE  computer  code  when  excited  by 
the  same  source.  This  chapter  presents  the  first  known  validation  of  a  time  dependent 
elastic  3D  finite  element  code  using  a  discretely  convolved  analytic  solution. 
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X 


Figure  4.  This  graphic  depicts  the  domain  of  a  semi  infinite  region  or  half¬ 
space.  The  x  —  y  plane  extends  from  — oo  to  +cx)  whereas  the  domain  of  0 
extends  from  the  origin  to  —  00  only. 

B.  ANALYTIC  FORMULATION 

Analytic  formulation  of  an  elastic  half-space  consists  of  deriving  the  elastic 
partial  differential  equation  governing  wave  motion,  determining  some  method  by 
which  to  solve  the  PDE,  and  finally,  finding  some  way  of  extending  the  solution  to  as 
general  of  a  forcing  function  as  possible. 

1.  Derivation  of  Elastic  Partial  Differential  Equations 

In  an  elastic  isotropic  medium  where  the  perturbing  influences  of  the  surface 
can  be  ignored,  waves  propagate  in  the  form  of  infinitesimal  stress  disturbances  over 
a  medium  [Ref.  20].  A11  idealized  analysis  ignores  the  static  body  force  of  gravity, 
and  assumes  that  the  entire  medium  has  a  uniform  density.  Figure  4  is  a  graphical 
depiction  of  this  phenomenon.  Of  interest  in  an  excited  elastic  material  are  displace¬ 
ment,  stress,  strain,  and  body  forces.  When  body  forces  are  absent,  the  development 
of  governing  equations  for  the  propagation  of  a  linear  elastic  wave  in  a  solid  material 
begins  with  the  equilibrium  equations  of  motion.  In  a  Cartesian  coordinate  system, 
following  the  convention  used  by  AchenbachfRef.  5]  and  Graff  [Ref.  4],  the  equilibrium 
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equation  of  motion  is  expressed  as 


82u 

PW 


Vr, 


which  translates  to 


f)2V  r)  7- 

V  LLX  U  I  XX  u  1  xy  u  1  xz 

P  dt2  dx  dy  dz 

f)2U  f)  J-  f)'!-  f)'T 

P  dt 2  dx  dy  dz 

d^u  <9t 

u  11  z  _  u  1  zx  u  1  zy  u  1  zz 

P  dt 2  dx  dy  dz 

and  can  be  written  compactly  in  tensor  form  as 

d2ut 

p  dt 2  TijJ ’ 


(in) 

(II-2) 
(II  .3) 

(HI) 

(II  .5) 


where  p  is  the  mass  density  of  the  material,  Ui  =  ( ux,uy,uz)T  is  the  displacement 
vector,  and  rtJ  is  the  stress  tensor.  Tensor  notation  is  a  concise  way  to  express  the 
continuum  mechanics  of  elastic  media  and  will  be  employed  throughout  the  remain¬ 
der  of  this  dissertation  unless  indicated  otherwise.  In  short,  a  single  index  denotes 
a  vector,  i.e.,  ut,  and  two  indices  will  denote  a  tensor  of  order  two  (a  matrix),  for 
example,  t13  .  Since  Cartesian  coordinates  are  used,  all  indices  denote  Cartesian  com¬ 
ponents  such  that  ui  =  ux,ii2  =  uy,U3  =  uz,ti2  =  rxy  etc.  A  summation  convention 
is  assumed  for  double  indices  so  that  whenever  a  subscript  appears  exactly  twice  in  a 
given  term,  that  subscript  will  take  on  the  values  1,2,3  successively,  and  the  resulting 
terms  summed.  The  summed  subscripts  are  simply  placeholders  or  dummy  variables 
since  it  is  immaterial  which  letter  is  used.  Note,  however,  in  clastic  theory  there  are 
terms  that  have  more  than  one  pair  of  dummy  indices  [Ref.  21],  Partial  differenti¬ 
ation  with  respect  to  a  Cartesian  coordinate  is  denoted  by  a  comma  preceding  the 
index  [Ref.  22],  Thus 


u 


hj 


duj 
dxj  ’ 


ij,k 


d'dij 

dxk 


(II.6) 
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From  these  rules  follow  that,  for  example,  ul>t  is  the  divergence  of  the  vector  rq  i.e., 
| +  +  |«3 1  while  Ui=  9n  [S  the  gradient  vector 

In  order  to  solve  11.1,  stress  must  be  expressed  in  terms  of  strain.  Hysteretic 
analysis  of  stress  and  strain  shows  that  a  one  to  one  relationship  exists  between  stress 
and  sufficiently  small  strains  in  an  elastic  body  [Ref.  23].  Therefore,  we  can  obtain 
linear  relations  between  stress  and  strain.  Equation  11.7  is  the  constitutive  equation 
representing  Hooke’s  law  [Ref.  24]  that  states  stress  is  proportional  to  strain 


T  Tjj  Cijkl^kl  4“  2/T6jj, 


(11.7) 


where  is  the  fourth  order  elasticity  tensor,  whose  elements  are  {i,j,  k,l  —  x,  y,  z ), 
€ij  is  the  strain  tensor,  A  and  y  are  modulus  constants,  and  Sij  is  the  well  known 
second  rank  tensor  known  as  the  Kronecker  delta,  whose  components  are  defined  as 


1 1  if  i  =  j 
1 0  if  i^j 


(II.8) 


We  can  relate  the  small-strain  tensor  to  displacements  within  the  limitations  of  the 
linearized  theory  of  deformation  by 


et,--_(utJ  +  uLk). 


(II.9) 


Substituting  the  strain-displacement  relations  into  Hooke’s  law  yields 


T{j  ^^k,k^ij  "b 


(11.10) 


which  is 


Tij  =  A(wi,i  +  u2, 2  +  u3i3)  +  y(uij  + 


(11.11) 


Because  of  the  symmetry  of  wave  propagation  in  an  isotropic,  homogeneous  media, 
ui,j  —  uj,i-  Therefore, 


Tij  —  A('Wii  +  u2,2  +  ^3,3)  +  2  yuij. 


(II-12) 
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It  follows  then  that  by  making  the  substitution  into  the  equation  of  motion,  and  using 
the  fact  that  the  elastic  stress  tensor  is  symmetric,  the  elastic  wave  equation  emerges 
as 

P~Q^  =  +  (A  +  H )ujdi .  (11.13) 

For  the  purpose  of  boundary  termination  the  frequency  domain  wave  equation  be¬ 
comes  useful  and  is  Fourier  transformed  to 

-pLU2Ui  =  pUijj  +  (A  +  n)ujS.  (11.14) 

so  that  in  a  homogeneous,  isotropic  medium,  this  equation  permits  plane-wave  solu¬ 
tions  of  the  form 

Um(xj,u;)  =  (11.15) 

where  A,m  represents  the  amplitude  and  the  polarization  of  the  plane  wave,  /q  is  its 
wave  number  in  cartesian  coordinates,  and  Xi  is  the  position  vector. 

2.  Method  of  Solution 

As  many  as  three  or  four  different  types  of  waves  may  exist  in  an  isotropic, 
homogeneous,  clastic  media  where  interfaces  and/or  free  surfaces  are  present.  This 
contributes  to  the  relative  complexity  of  elastic  solid  wave  problems  compared  to 
equivalent  problems  in  acoustics  and  electromagnetics.  Where  no  body  forces  exist, 
11.13  is  the  displacement  equation  of  motion  and  represents  a  system  that  has  the 
disadvantageous  feature  of  coupling  the  three  displacement  components  together.  To 
solve  this  directly  would  require  solving  a  sixth  order  partial  differential  equation.  A 
more  convenient  and  commonly  used  method  for  solving  the  equation  would  be  to 
express  the  components  of  displacement  in  terms  of  their  potentials  and  derivatives 
of  potentials  [Ref.  5].  Consider  a  decomposition  of  the  displacement  vector  of  the 
form 

Ui  T  Gipq'lpq  p  (11.16) 
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where  0  is  referred  to  as  the  scalar  potential  or  irrotational  portion,  curl  ip q  is  referred 
to  as  the  vector  potential  or  rotational  portion,  and  eipq  is  a  common  and  frequently 
used  rank  three  alternating  tensor,  whose  components  are  as  follows: 


+1  if  ipq  represents  an  even  permutation  of  123 


e.  = 
z'ipq  s 

0  if  any  of  the  indices  are  equal 

—  1  if  ipq  represents  an  odd  permutation  of  123 

(H-17) 

By  this  decomposition  primary  variables  are 

Ul  =  0,1  +  03,2  -  02,3 

(11.18) 

U2  =  0,2  +  01,3  -  03,1 

(11.19) 

U3  =  0,3  +  02,1  -  01,2 

(11.20) 

The  advantage,  of  course,  of  this  convention  is  that  it  allows  us  to  separate  the  dis¬ 
placement  field  into  two  distinct  components,  compressional  waves  and  shear  waves. 
We  can  use  this  decomposition  and  solve  each  potential  individually  then  reconstruct 
to  solve  the  entire  system.  We  start  by  substituting  11.16  into  11.13  to  yield 


P(0,i  T  &ipq'ipq,p)  1^(4, i  T  ^ipq^Pq,p)  ,jj  T  (A  T  p)(0j  T  ^ jpq'P q.p)  ,ji-  (H-21) 

which  is  further  evaluated  as  (see  [Ref.  21]) 

P(0,i  T  &ipq'ipq,p)  P0, ijj  T  P^ipq'^Pq:pji  T  (A  T  /00,p'j  T  (A  T  l^)&jpq1pq,pji-  (11.22) 

Since  eTpqipqmi  =  0  we  obtain 

P(0,i  T  f'ipq'ipq,p)  P ( 4.ijj  T  ^ipq'*Pq;pjj^)  T  (A  T  p)  0,pj  (11.23) 

and  upon  rearranging  terms  we  have 

((A  T  2/i) (pjj  P0),i  T  ^"ipq(l^4q,jj  P0<?),P  0.  (11.24) 

We  can  see  now  that  11.16  satisfies  the  equation  of  motion  if 

(A  +  2 n)(ptjj  -  p4>  =  0  (11.25) 
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Compressional 
or  P-wave 


Shear  or 
S-wave 


Figure  5.  This  graphic  is  a  common  depiction  used  to  represent  the  propa¬ 
gation  of  P-waves  and  S-waves  in  an  elastic  body.  The  behavior  of  P-waves 
involve  compressional  or  dilatational  motion  that  changes  the  volume  of 
the  material  as  the  wavefront  passes  while  S-wave  motion  represents  a 
shearing  of  the  material,  but  does  not  change  the  material’s  volume  as  the 
wave  progresses. 


and 


/# q,jj  ~  Pi>q  =  0. 


(11.26) 


11.25  and  11.26  are  the  uncoupled  wave  equations  of  an  elastic  media  and  reveal  the 
total  motion  as  being  composed  of  a  curl-free  pressure  wave  (P-wave)  traveling  with  a 
speed  2 ,  and  a  divergence-free  shear  wave  (S-wave)  traveling  with  a  speed 2 . 

By  denoting  the  wave  speeds  as  cp  and  cs  respectively,  the  equations  of  motion  become 


c2p(t){xi,t)dj  =  4>{xi,t )  (11.27) 

and 

C2s^{xh  t)qjj  =  V>(£i,  t)q.  (11.28) 

Equation  11.27  models  the  behavior  of  P-waves  which  involve  compressional  or  di¬ 
latational  motion  that  changes  the  volume  of  the  material  as  the  wavefront  passes 
while  equation  11.28  models  S-wave  motion  representing  a  shearing  wave  that  does 
not  change  the  volume  of  the  material  as  the  wave  progresses.  Figure  5  is  a  common 
depiction  used  to  graphically  represent  these  two  motions  on  an  clastic  body.  If  we 
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take  proposed  decomposition  11.16  and  place  it  back  into  II.  10,  the  stress  components 
can  be  written  in  terms  of  their  displacement  potentials  as 

7~ij  A[0,fcfc  “F  (^'kpq^Pq/p)  ,0  Aj  ~F  p  [0,jj  “0  (c  ipq^fiq,p)  J  4*,ji  "I-  (ejpq'tP<hp) >*]  (H-29) 

where  symmetry  yields 

T~ij  A[0,/jfc  "F  (Cfcpq0</,p)  ,k\ &ij  "F  2  p  [0,jj  ~F  (Cipq0q,p)  j]  •  (11.30) 

For  example,  suppose  we  wanted  to  know  the  stress  component  Tn,  then 

Tn  =  A[0,n  +  (eip?0g,p),i  +  0,  22  +  {e2pq'lPq,  p),  2  +  0,33  +  (esp^^p)^]  An 
~F  2/1  [0,11  ~F  (Clpg0g,p)  ,l]  • 

After  summing  over  the  p  and  g  we  have 

Til  =  A  [0,11  +  (ei2303,2),l  +  (^  13202,3  )  ,l]  All  +  A[0,22  +  (e21303,l),2  +  (^23101, 3)  ,2]  All 
+  A [0,33  +  (^32101, 2), 3  +  (e31202,l),3]^ll  +  2 p  [0,n  +  (ei2303,2),l  +  (^13202, 3) ,l]  • 

By  taking  into  account  the  even  and  odd  permutations  of  eipq  according  to  equation 
11.17  and  Kronecker’s  delta ,  we  have 

Tu  =  A  [0.11  +  03,21  —  02, 3l]  +  A  [0,22  —  03,12  +  01,32]  +  A  [0,33  —  01,23  +  02,13] 

+  2 p  [0,11  +  (03,2  —  02, 3), l] 

Tn  =  A[0.n  +  0,22  +  0,33  +  (03.1  _  03, 1)2  +  (02,3  _  02, 3)1  +  (01,2  —  01, 2)3] 

+  2/i  [0,n  +  (03,2  -  02, 3), l]  • 

Judicious  collection  of  the  0  terms  reveals  the  cancelation  of  terms  multiplied  by  A, 
simplifying  to 

Tn  =  A0,fcfc  +  2 /I  0,ii  +  (03,2  —  02,3)  !  , 
and  accordingly  for  the  remaining  stress  components  we  have, 

T22  =  A 0^^,  +  2/i  0,22  —  (03,i  —  0i, 3), 2 
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H(t) 


-z 

Figure  6.  Point  load  problem  in  Cartesian  coordinates. 

733  =  A0,fcfc  +  2 H  0,33  +  ('02,1  —  01,2)  3 

T12  =  T21  =  p  [20,12  +  (03,2  -  02,3) ,2  ~  (03,1  “  01,3) ,1 

7-23  =  7-32  =  Id  [20,23  “  (03,1  “  01, 3^3  +  (02,1  ~  01, 2)^ 

7-31  =  T13  =  [I  20,31  +  (03,2  -  02, 3), 3  +  (02,1  -  01,2),!  • 

Thus  by  solving  11.27  and  11.28  respectively,  we  can  find  not  only  displacement  from 
11.16,  but  also  the  amount  of  stress  at  any  point  as  well. 

a.  The  Half- Space 


In  order  to  formulate  a  problem  for  an  elastic  half-space  (z  <  0),  we 
need  boundary  conditions  at  the  surface  (see  Figure  6).  Our  source  will  be  a  point 
load  on  the  surface  z  —  0,  and  boundary  conditions  will  be  expressed  as  components 
of  stress.  Thus, 

TijUj  =  0  (11.31) 


except 


733  =  -WH{t)8{x)8{y) 


(11.32) 


where  H(t)  is  the  time  dependent  Heaviside  function  and  W  is  a  force.  Together,  with 
the  spatial  two  dimensional  delta  function  S(x)S(y),  r33  specifies  a  pressure  or  traction 
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on  the  surface.  The  problem  is  not  completely  formulated,  however,  until  we  make  a 
statement  about  the  initial  conditions.  In  this  case,  our  medium  is  undisturbed  until 
some  time  t  where  it  is  excited  by  a  pressure  directed  downward.  This  gives  initial 
conditions  as 

(j)(xi,  0)  =  tj>(xi,  0)  =  0  (11.33) 

^k(xi,  0)  =  ipk(xi,  0)  =  0  (11.34) 

Taken  all  together  11.33,  11.34,  along  with  11.27  and  11.28  form  a  well  posed  system 
of  PDE’s  which  can  be  solved  individually  to  find  all  components  of  displacement  in 
the  elastic  media. 

b.  Integral  Transforms 


The  goal  of  this  section  is  to  obtain  integral  representations  of  the  po¬ 
tentials.  These  potentials  can  then  be  used  in  conjunction  with  the  stress  representa¬ 
tions  to  find  unique  solutions  to  the  potentials  and  thereby  determine  the  half-space 
surface  displacements.  Using  the  one-sided  Laplace  transformations  [Ref.  25],  where 
Br  denotes  the  Bromwich  inversion  path  in  the  complex  p-plane, 


—  /*00 

f(r,z,u)=  /(r,z,f)e“^dt,  (11.35) 

J  o 

f(r ,z,t)  =  — —  f  J(r ,  z,u>)eu>tdw,  (11.36) 

ZlTl  J  Br 

all  time  dependence  of  equation  11.27  and  equation  11.28  are  removed.  Thus, 

u2  - 

(f>(xi,u)tu  =  —4>{xi,u)  (11.37) 

and 

uj2  - 

i/)(xi,u;)k,ii  =  -~il)(xi,L;)k  (11.38) 

u s 

are  the  integral  transform  representations  for  the  two  potential  equations.  The  next 
step  uses  a  Double  Fourier  Transform  pair  [Ref.  26],  given  as 

~  1  r°°  r00 

TX6,6)  =  ^/  /  f(x1,x2)e*lXle*aX*dx1dx2  (11.39) 

Z7T  J— oo  J— oo 
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f(x,y) 


(11.40) 


1  POO  POO  ~ 

=  7T  / 

Z7T  J  —  oo  J  —  oo 

over  the  entire  infinite  xi(x)  space  and  x2 (y)  space.  Note,  also,  that 

e*£i xei&y  _  ei{t,ix+&y)  _  e*€i-a;i 


(11.41) 


The  x3{z)  dimension  remains  untransformed  in  the  analysis.  This  converts  the  PDEs 
into  ODEs  with  respect  to  x 3.  For  <j)(xi,uj)  and  ^(ay,  to)k  we  have: 

u2  - 

0(xi,  x2,  x3,  to)  a  =  -ir<t>(xi,  x2,  x3,  to)  (11.42) 

Lp 

u2  - 

il>(xi,  x2,  X3,  u)ktu  =  -fil>(xi,  x2,x3,  u)k.  (11.43) 

After  application  of  spatial  transforms  in  the  x\{x),  x2 (y)  coordinates,  and  integration 
by  parts  we  derive 


0(£i,  £2,  %3,  <^),33  —  +  (£1  +  £2)^  0(£i>  £2,  ^3,  u)  (11.44) 

^(£1,  £2,  x3,  u)kt 33  =  +  (£1  +  £2)  j  ^(£1,  £2, 2:3,  u)k.  (11.45) 

The  above  equations  are  now  second  order  ODEs  with  respect  to  £3(2:).  Solutions  of 
equation  11.44  and  equation  11.45  take  the  form 


0(£l,  £2,  X3,  to) 

^(£l,  £2,2:3 ,w)k 


*h(£i,  £2,  <^’)e± 

^(£i,£2,w)fce 


a +<«+«?))’ 

±(f+<e;+®) 


013 


1 

013 


(11.46) 

(11.47) 


At  this  point  it  is  interesting  to  note  that  in  a  homogenous  isotropic  material  clastic 
waves  propagate  equally  well  in  all  directions.  This,  of  course,  is  intuitive  since  the 
material  acts  on  the  speed  of  the  wave  in  the  same  way  at  every  point  as  it  propagates 
through  the  substance.  Therefore,  the  motion  is  invariant  with  respect  to  £2  if  the 
wave  is  traveling  normal  to  £1  and  vice  versa.  This  means  that  the  complexity  of  the 
problem  can  be  reduced  significantly  by  using  cylindrical  polar  coordinates  (r,z,9). 
Figure  7  shows  the  mapping  of  polar  (i.e.,  (2  =  +  £|,  z  =  x3)  instead  of  Cartesian 
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H(t) 


-z 


Figure  7.  Point  load  problem  in  polar  coordinates. 


coordinates  where  the  subscript  r  in  £r  denotes  radial  motion  in  the  x  —  y  plane.  For 
convenience,  £  =  £r.  Primary  variable  equations  11.18,  11.19,  and  11.20  are  converted 
to  their  polar  counterparts  as 


Ur  <j),r  ^z,9  ^Pd,z 

Ue  =  4>,e  + 

Uz  =  4>,z  +  ^9,r  -  4>r,6 


1 

"Si  ^ 

+ 

~0- 

(11.48) 

1  d(j)  dTpr  dl^z 
r  dO  dz  d9 

(11.49) 

d(p  |  1  dfyor)  1  ch/y 
dz  r  dr  r  dO 

(11.50) 

Since  all  9  dependence  can  be  in  effect  removed  because  wave  motion  is  axially  sym¬ 
metric.;  see  [Ref.  27]  and  [Ref.  5],  we  are  essentially  left  with 


Ur  —  0,r  ^p9,z 

Uz  (fr,z  "b  IpO, r 


dtp  d^Q 
dr  dz 
d<j)  1  d{ij)  or) 

dz  r  dr 


(11.51) 

(11.52) 


Observe  that  now  we  are  actually  manipulating  two  uncoupled  problems  of  wave 
motion.  ur  and  uz  depend  only  on  the  quantities  (j)  and  ipo  which  are  governed  by  the 
independent  scalar  wave  equations  11.53  and  11.54.  We  say  equation  11.54  is  scalar 
because  horizontal  and  vertical  displacement  equations  are  computed  from  only  the 


20 


ipe  component  of  the  vector.  Thus, 


U! 


(/>(£,  Z,v),zz=  +  £  \(f>{i,z,u) 


UJ 


^(^,z,uj)e,zz=  +  £  m£> 


becomes  the  two  ordinary  differential  equations  of  potentials 


(11.53) 

(11.54) 


where  a  and  (3  are 
will  be  of  the  form 


d2cf) 
dz 2 


(11.55) 


d?j>e 

dz 2 


=  f32^e 


(11.56) 


+  £2  and  +  £2  respectively.  Solutions  to  11.55  and  11.56 


^z,u)  =  <$>{t,u)e±az  (11.57) 

^z,u)  =  ^u)e±fiz  (11.58) 


With  respect  to  the  new  system  of  cylindrical  polar  coordinates  ( r,0,z ),  the  axially 
symmetric  motions  in  a  half-space  break  down  to  the  following  stress  equations 


tzz  =  (A  +  2/i)uZiZ  +  A  r  1(rur)jr  (11.59) 

Tzr  di.^r,z  T  UZlr)  (11.60) 

where  again  we  have  (rzz,  rzr )  as  the  components  of  the  stress  tensor.  Next  the  bound¬ 
ary  conditions  on  the  surface  have  to  be  converted  to  cylindrical  polar  coordinates. 
Thus, 

rzz  =  -WH{t) ^  (11.61) 

and 

Trz  =  0  (11.62) 

are  converted  to  their  polar  equivalents.  Initial  conditions  properly  converted  simply 
become 

0(r,  z,  0)  =  0(r,  z,  0)  =  0  (11.63) 
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ij}(r,  z,  0)  =  ^(r,  z,  0)  =  0 


(11.64) 


Now  we  use  a  couple  of  transform  pairs  to  seal  the  deal.  For  time, 

we  keep  the  one 

sided  Laplace  transform  11.35  and  11.36.  In  place  of  the  double  Fourier  transform 

(11.39  and  11.40)  used  earlier,  the  Hankel  transform  is  enlisted  to  do  the  task  because 

of  its  axisymmetry  and  use  of  Bessel  functions.  The  appropriate  definitions  are  as 

follows  [Ref.  25] 

_  roo 

L{f}  =  f(r,  z,u)  =  /  f(r,  z,  t)e~utdt, 

Jo 

(11.65) 

L{f}  =  /(A  z,t)  =  J-  [  /(A  A  u)(?*dwt 

Zj[X  J  Br 

(11.66) 

_  _  rOO  _ 

Hu{f}  =  f(£,z,u)  =  /  f(r ,  z,u))Ju(£r)rdr , 

Jo 

(11.67) 

~  _  roo  _ 

Hv{f}  =  f(r,z,u)  =  /  f(£,z,u)J„(€r)€d£ 

Jo 

(11.68) 

where  equation  11.65  and  equation  11.67  provide  the  direct  transforms  and  equation 
11.66  and  equation  11.68  yield  inversions.  Br  denotes  the  Bromwich  inversion  path 
in  the  complex  p-plane,  and  Jv()  are  Bessel  functions  of  the  first  kind  of  order  v. 
Applying  the  transforms  to  the  displacements  yield 


r  oo 

v,r(€,z,uj)=  /  ur(r,  z,cj)J1(£r)rdr 

’) = r  -  it) M(r)rdr 

w)  =  jf  Ji(£r )rdr  -  jf  "  j  -J\  (£r)rdr 

r  oo  ^  /'oo 

=  -f  /  <f>Jo(£r)rdr  -  —  /  ^eJ\{ir)rdr 
Jo  dz  Jo 

ul(£,z,w)  =  -£4>°(£,z,u)  -  —fr(£,z,u) 
for  radial  displacement  and 

roo 

u°z(£,,z,uj)  =  j  uz(r,  z,u)Jo(£r)rdr 

^z-°j)=r{fz+^+T)M(r)rdr 


(11.69) 

(11.70) 

(11.71) 

(11.72) 

(11.73) 

(11.74) 

(11.75) 
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«°(f,  z,u)  =  j  (jjz^  Ja{£r)rdr  +  ^  J0(£r)rdr  (11.76) 


d  r°°  /*°° 

u°g(£,z,u)  =  —  y  <j)J0(£r)rdr  +  £  jf  i/jeJi(^r)rdr 
«°(f ,  z,  u)  =  ^0°(£,  z,  v)  +  ,  2,  w) 


(11.77) 

(11.78) 


for  vertical  displacement.  The  necessary  stresses,  equations  11.59  and  11.60,  are  re¬ 
placed  with  the  equivalent  potentials  and  then  similarly  transformed 


TZz  (A  T  2/0  / 0  zz  T  0  2r  T  0,2^  T  (0  r  0  2:  T  T‘(p ,rr  ^"0zr) 

L  r  J  r 

Tzz  =  (A  +  2/x)0°2z  +  2/ii/o{0,r},2  +  2/t.f/o  j“0,2  j  —  £"A0° 


from  equation  11.55,0°^  =  a20° 


=  (A  +  2/0C000  +  2/x77o{0,r}>2  +  2/itf0  ^  -0,2  )>  -  A0° 


^22  —  (A  +  2 /x)  +  02J  0°  +  2//77o{0,r},2  +  2 0^0,2  —  £2A0 

<90 1 


2\I0 


^zz  ~  d  —  +  2£z  )  0U  +  2£/t- 


cn 


<9;£ 


and 


Tzr  =  -2^^-  -  /*  01 


<9z 


Now  we  transform  the  boundary  conditions  on  the  surface  to  become 

W 


T°  = 


2ttuj 


and 


rl  =  0. 


(11.79) 

(11.80) 

(11.81) 

(11.82) 

(11.83) 

(11.84) 

(11.85) 

(11.86) 


Equations  11.57  and  11.58  are  assumed  to  be  the  solutions  to  0°  and  01.  Therefore, 
placing  the  assumed  solutions  into  equations  11.83  and  11.84  with  boundary  conditions 
11.85  and  11.86  yields  a  system  of  stress  equations  which  can  be  used  to  determine  <f> 
and  T.  Thus, 

I*  +  2<0  *(£,w)e-“  +  2^0  (11.87) 
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and 


-2£/i£  ($(£,  uj)e-az)  -  /i  +  2Z2]j  “)e~Pz  =  0  (H.88) 


which  when  evaluated  at  the  surface  z  —  0  gives  the  two  equations 

W 

2ttlu 


ui 


u> 


R  -^  +  2^  $-2£//  — +  T  4>  = 


and 


The  solutions  are 


-2f  w  +  r  $  +  -v  +  2r  *  =  o- 


*=-m 

V  27 r/icu  / 


2  , 


2?  + 


2^  +  s  -4 eJe  +  (£Je  +  i4 


and 


^  f  IT  ^  _ 2^C2  +  % _ 

^+iy_4^v/{2+|v^T 


(11.89) 


(11.90) 


(11.91) 


(11.92) 


This  now  allows  the  determination  of  explicit  vertical  and  radial  displacements  ([Ref. 
5])  by  using  equations  11.73  and  11.78,  such  that 


and 


uj:(€,z,w)  =  -f$(f,u>)e  az  -  ^(^,u)e  9z 


u%,z,u)  =  j-^{^)e~az  +  ^{^)e-9z. 


(11.93) 


(11.94) 


Analytic  Solution 


Integral  inversions  of  equations  11.93  and  11.94  according  to  integral  in¬ 
version  equation  11.68  produce  the  following  equations  (see  [Ref.  27])  for  the  Laplace- 
transformed  displacements  at  the  surface,  i.e.,  z  —  0 


ur(r,  0,u>)  = 


W  r°°  ¥  (2V^  +kp\Jt2  +  ks  ~ks-  2^2)  Ji(Zr) 


2nyu  Jo 


R(£,u) 


(11.95) 
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Pekeris  Closed  Form  Solution 


Figure  8.  Vertical  surface  displacement  uz(r,  z,  t )  due  to  a  concentrated  point 
source  directed  downward  on  the  surface.  Poisson’s  ratio,  is,  is  |  with  A  =  /i. 
The  measurement  is  located  one  unit  radially  from  the  point  of  impact.  P 
denotes  the  arrival  time  of  compressional  waves,  and  S  of  the  shear  wave. 
The  Rayleigh  wave  propagates  at  the  speed  of  the  singularity. 


uz(r,0,u) 


W  (sjt2  +  kp )  %J0(£r) 

2,71  /jluj  Jo  R(£,  uS) 


where  R(£,u>)  is  the  Rayleigh  function  defined  as 


r((,v)  =  <2f2 + kif  -  4 


(11.96) 


(11.97) 


and  variables  kp,  ks  are 

kp  =  —  (11.98) 

Cp 

and 

k8  =  —  (11.99) 

Cs 

respectively.  The  challenge  is  to  perform  the  transform  inversion  of  the  displacement 
integrals.  Many  have  accomplished  this  with  various  techniques.  There  have  been 
numerical  approaches  [Ref.  28]  as  well  as  closed-form  solutions.  Most  notable  of  the 
closed-form  solutions  would  be  that  of  Pekeris  [Ref.  2],  Achenbach  [Ref.  5]  con- 
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Pekeris  Closed  Form  Solution 


Figure  9.  Horizontal  surface  displacement  ur(r,z,t)  due  to  a  concentrated 
point  source  directed  downward  on  the  surface.  Poisson’s  ratio,  v,  is  | 
with  A  =  /i.  The  measurement  is  located  one  unit  radially  from  the  point 
of  impact.  P  denotes  the  arrival  time  of  compressional  waves,  and  S  of  the 
shear  wave.  The  Rayleigh  wave  propagates  at  the  speed  of  the  singularity. 

sidered  the  anti-plane  shear  version  of  the  problem  and  employed  cylindrical  integral 
methods  which  take  advantage  of  the  wave  propagation’s  symmetry.  Both  Pekeris 
and  Achenbach  used  the  Cagniard-de-Hoop  method  as  an  analytical  technique  for 
evaluating  the  multidimensional  Fourier  integrals,  and  obtained  an  inversion  proce¬ 
dure  for  the  displacement  integrals  which  then  could  be  solved  using  partial  fraction 
decomposition.  If  the  radial  distance  from  the  point  source  is  unity  and  Poisson’s  ra¬ 
tio  is  |  with  A  —  p,  then  the  vertical  and  horizontal  displacement  of  a  particle  at  the 
surface  is  displayed  in  Figures  8  and  9.  Figure  8  depicts  the  vertical  surface  displace¬ 
ment  uz(r,z,t)  due  to  a  concentrated  point  source  directed  downward  on  the  surface 
while  Figure  9  shows  the  corresponding  horizontal  surface  displacement  ur(r,  z,t ).  In 
both  figures  the  longitudinal  wave  arrival  time  is  clearly  seen.  The  shear  wave  ar¬ 
rival  (though  not  as  obvious)  is  also  visible  by  a  sudden  change  in  the  velocity  of  the 
wavefront  as  it  passes  the  observation  point.  The  Rayleigh  wave  arrival  time  is  purely 
theoretical  and  coincides  with  the  speed  of  a  singularity  that  propagates  along  the  free 
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surface.  It  is  also  noted  that  because  of  the  functional  singularity,  the  amplitude  of 
the  Rayleigh  wave  cannot  be  accurately  determined.  Technically,  the  singular  point 
has  an  amplitude  of  infinity.  Certainly,  in  practice  the  amplitude  is  finite.  What 
Pekeris’  closed-form  solutions  provides  are  the  arrival  times  at  which  all  three  wave 
fronts  pass  an  observation  point.  In  this  case  that  point  is  1  meter  radially  from  the 
source.  This  is  extremely  useful  in  helping  to  determine  if  the  numerical  FE  model 
has  the  proper  phase  speeds  within  the  computational  area.  So  although  the  ampli¬ 
tude  of  the  surface  wave  may  not  be  easily  expressed  in  terms  of  a  definite  magnitude, 
the  arrival  time  of  a  disturbance  propagating  in  the  domain  can  be  a  key  indicator 
of  its  presence  in  a  FE  model.  The  distribution  of  total  energy  for  a  single-element 
radiator  among  the  shear,  compressional,  and  Rayleigh  waves  are  25.8  percent  for 
shear,  6.9  percent  for  compressional,  and  67.4  percent  for  Rayleigh  respectively  as 
computed  by  Miller  and  Pnrsey(1955)  [Ref.  29]  for  the  idealized  Poisson’s  ratio  | 
with  A  =  /i,  and  experimentally  observed  by  Wood(1968)  [Ref.  30].  Figure  10  is  a 
slice  of  SAFE-T’s  FE  model  which  graphically  shows  the  different  wave  propagation 
speeds  as  well  as  different  displacement  amplitudes  among  the  three  waves  present 
in  a  seismic  event.  The  anisotropic  propagation  in  the  —  z  direction  is  because  the 
elements  are  not  square.  This  means  that  the  shape  of  the  elements  must  be  ”iso”  if 
one  wishes  to  model  isotropic  behavior. 

3.  Extension  to  Arbitrary  Transient  Source  Wave  Form 

Transient  sources  such  as  the  Heaviside  function  or  the  Dirac  Delta  function 
allow  for  the  analytic  solution  as  given  above  for  Lamb’s  problem,  but  they  are  not 
realistic  because  of  the  non-physical  singular  impulse  associated  with  them.  How¬ 
ever,  the  Heaviside  function  as  a  source  and  analytically  computed  by  Pekeris  can 
be  extended  by  the  principle  of  superposition  and  integral  convolution  to  model  any 
arbitrary  transient  waveform  as  an  impact  source.  The  necessity  of  such  a  convention 
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Rayleigh  wave 


Figure  10.  This  slice  of  SAFE-T’s  FE  model  graphically  shows  the  differ¬ 
ent  wave  propagation  speeds  as  well  as  different  displacement  amplitudes 
among  the  three  waves  present  in  a  seismic  event.  The  anisotropic  prop¬ 
agation  in  the  —  z  direction  is  because  the  elements  are  not  square. 

is  seen  in  the  fact  that  it  is  impossible  to  numerically  model  the  instantaneous  change 
of  state  associated  with  the  Heaviside  and  Dirac  Delta  functions. 

a.  The  Source 

Physically  speaking,  sources  can  be  similar  to  near-perfect  adhesions. 
One  example  of  this  would  be  a  weight  dropped  onto  damp  soil.  Sources  can  also 
be  simulated  as  near-perfect  rebounds  like  a  steel  pellet  onto  some  type  of  granite 
surface.  Discussions  on  this  matter  with  solid  mechanical  engineers  along  with  the 
work  of  Rumph  [Ref.  9]  has  led  the  author  to  select  a  single  waveform  which  will  be 
suitable  for  all  models  in  this  report.  The  source  used  in  Rumph’s  experiment  is  a 
Haversine  pulse.  The  Haversine  waveform  is  produced  by  shifting  the  phase  of  the 
sine  wave  by  90  degrees  and  then  adjusting  the  offset  to  set  the  baseline  to  zero.  A 
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f  (t) 


Figure  11.  Gaussian  surface  point  source  is  modeled  after  Rumph’s  experi¬ 
mental  Haversine  source.  The  Haversine  waveform  is  produced  by  shifting 
the  phase  of  the  sine  wave  by  90  degrees  and  then  adjusting  the  offset  to 
set  the  baseline  to  zero.  This  Gaussian  has  similar  properties  as  the  Haver¬ 
sine  waveform  source,  is  easier  to  implement  into  numerical  code,  and  is 
everywhere  differentiable. 


f  (t) 


Figure  12.  Derivative  of  Gaussian  surface  source  shows  that  the  derivative 
is  continuous  everywhere. 
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simpler  source  to  use  is  a  Gaussian  source  shown  in  figure  11  modeled  with  the  same 
characteristics  as  Rumph’s  experimental  Haversine  pulse.  It  starts  with  zero  slope, 
passes  smoothly  through  a  peak  value  and  returns  symmetrically  to  zero  again.  Since 
it  is  a  Gaussian,  it  has  a  first  derivative  which  is  continuous  everywhere,  and  allows 
for  easy  modifications  to  its  amplitude,  width,  and  location  of  peak  (see  figure  12) 
accomplished  through  the  function 

/ t—peak \ ^ 

f(t)  —  (Amplitude) e  v  width  >  (II. 100) 


with  derivative 


df(t) 

d.t 


-2 


Amplitude  \ 

(width)2  J 


e 


t—peak \ ^ 

width  >  (t  —  peak ) 


(II. 101) 


In  addition,  equations  II.  100  and  II.  101  allows  the  source  to  be  tailored  to  specific 
characteristics.  Namely,  by  adjusting  the  pulse  width,  one  can  change  the  power 
spectrum  of  the  pulse  generated.  The  power  spectrum  of  the  pulse  in  figure  12  is 
instrumental  in  determining  the  spatial  resolution  of  the  numerical  grid  which  in 
turn,  through  stability  requirements,  enforces  a  limit  to  the  time  step  that  can  be 
taken  between  each  calculated  solution. 

In  the  present  instance,  the  peak  power  of  the  pulse  occurs  at  599.675 
hertz  (see  figure  13)  which  corresponds  to  a  Rayleigh  wavelength  of  one  meter  given 
parameters  for  a  sandy  medium  (see  Appendix  C).  The  spatial  resolution  as  mea¬ 
sured  by  the  largest  distance  between  nodes  in  the  numerical  domain  is  less  than  A 
of  a  meter,  guaranteeing  sufficient  frequency  for  a  Rayleigh  wave.  As  there  are  other 
waves  propagating  in  the  elastic  medium,  they  should  also  be  weighed  in  mesh  res¬ 
olution  considerations.  Because  the  phase  speeds  for  compressional  waves  and  shear 
waves  are  greater  than  the  Rayleigh  wave  speed,  their  corresponding  wavelengths  at 
the  dominant  frequency  will  be  even  greater  than  those  of  the  Rayleigh  wavelength 
guaranteing  even  better  resolution. 
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F(OJ) 


Figure  13.  The  power  spectrum  of  the  derivative  of  Gaussian  surface  source 
reveals  the  dominant  frequency. 

b.  Convolution 


This  section  outlines  the  development  of  a  more  useful  analytic  result. 
Starting  with  the  specific  problem  treated  by  Pekeris,  expressions  for  vertical  and 
horizontal  displacements  appear  below  where  r  represents  the  dimensionless  quantity 
— ,  and  parameters  k  and  k  are  related  as 


k2  = 


(11.102) 


The  vertical  displacement  is  expressed  as 


V(r) 


fo 


\J  3\/3+5  ,  \J  3\/3— 5 

_ ZL  / R  \/ 3p3+5  ^ 

48  l 

_  7T 

8 


T  K  V3 
7s  <T<1 


1  <  r  <  7 


7  <  r 
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Applied  Force  Conversion  Table 

Component 

Unit  Step 

Impulse 

Arbitrary 

H(t) 

<5(r) 

/(T) 

Displacement  [uz  or  ur) 

V{t) 

dV(r) 

dr 

V(t) 

Velocity  {uz  or  ur) 

dV(r ) 
dr 

d2V{  t) 
dr2 

dV(r) 

dr 

Acceleration  (ilz  or  ur) 

c12V(t) 
cLt 2 

d3V(r) 
dr 3 

d2V(r) 

dr2 

Table  I.  Conversion  formulas  for  various  source  wave  forms  and  measured 
quantities  allow  the  use  of  integral  convolution  to  find  solutions  for  arbi¬ 
trary  transient  source.  Note:  (1)  V(r)  =  *V(r)  and  (2)  For  horizontal 

motion,  replace  V(t)  with  H(t). 


and  horizontal  displacement  as 


r  < 


d3 


H(t)  = 


si  {6*(*0  -  isn(sfc2,  k)  +  +  nl7V#V'}  7  <  r  <  1 

(20-12V3 
(6-4^3)- 


(6— 4v/3)-1  1  (6+4V3)- 

{6K(k)  -  1811(8,  k)  +  nl;”-^!!’K|  +  gjggwjM} 


l6\/i6 

Preceding+  yfr1 


A<r<1 


7 


7  <  T 


respectively  with  7  =  «  1.08766.  and  II (n,  k)  are  elliptic  integrals  of 

the  first  and  third  type.  By  convolving  an  arbitrary  transient  source  function  /(f) 
and  the  functions  given  above  using  linear  superposition  over  time,  it  is  possible  to 
find  accurate  comparisons  for  the  FE  model. 

Since  the  time  history  of  /(f)  is  known,  Mooney’s  [Ref.  6]  technique 
can  be  followed  and  the  use  of  the  Duhamel  integral  with  the  quantities  outlined  in 
Table  I  are  employed.  Digital  sampling  of  the  source  fit)  and  discrete  convolution 
are  used  to  calculate  the  response  of  the  free  surface  and  surrounding  media  due  to 
the  source  pulse  by  the  relation 


V(t) 


df  (r) 

dr 


*V(r) 


(11.103) 
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Surface  Response  1m  away 


time  (seconds) 


Figure  14.  Triple  view  of  the  two  waveforms  (2nd  and  3rd  graphs)  which 
combine  to  demonstrate  the  response  on  the  surface  1  meter  away  from 
the  source  (top  graph). 
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Surface  Response  1m  away 


Figure  15.  Solution  of  derivative  of  a  Gaussian  source  by  linear  superposi¬ 
tion  over  time  to  determine  the  surface  response  1  meter  away  from  the 
source  of  excitation. 

or  specifically, 

V(t)  =  y  ^>V(t  -  n).  (11.104) 

The  relationship  between  discrete  and  continuous  convolution  is  well  documented 
[Ref.  31].  In  fact,  the  most  important  applications  of  the  discrete  convolution  occur 
not  by  sampling  periodic  functions  but  rather  by  approximating  continuous  convolu¬ 
tions  of  waveforms  [Ref.  32],  Figures  14  and  15  display  convolutions  made  using  the 
derivative  of  a  Gaussian  for  fit)  of  the  form  of  equation  II.  100. 

C.  NUMERICAL  APPROACH 

1.  FE  Modeling  of  Partial  Differential  Equations 

A  benchmark  ’’analytic”  solution  to  the  well-posed  partial  differential  equa¬ 
tions  and  boundary  conditions  for  a  point/line  loaded  elastic  half-space  has  been 
discussed.  Except  for  a  few  simple  cases  it  is  nearly  impossible  to  find  analytic  so- 
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SAFE-T  Vertical  vs.  Exact  Results 
For  Unit  Step  Function 


Figure  16.  Solid  Adaptive  Finite  Element  -  Transient  (SAFE-T)  is  used  to 
model  vertical  displacement  due  to  an  instantaneous  Heaviside  change  of 
state.  92  percent  accuracy  is  achieved  with  2,764,800  degrees  of  freedom. 
At  time  steps  are  taken  at  0.0005  seconds.  This  particular  calculation  was 
done  on  a  Compaq  laptop  computer  with  an  AMD-64  Athlon  processor 
and  2  gigabytes  of  memory.  It  took  17.78  hours  of  wall  time  with  3.05 
hours  of  that  being  actual  CPU  time. 

lutions  to  the  governing  partial  differential  equations  with  nonsymmetric  geometries 
and  complex  boundaries.  The  finite  element  method  (FEM)  offers  an  alternative  to 
the  somewhat  limited  class  of  problems  for  which  analytic  solutions  can  be  found  [Ref. 
7].  By  replacing  the  differential  equation  with  an  equivalent  algebraic  system  of  equa¬ 
tions  i.e. ,  Au  =  /,  it  is  possible  to  produce  an  approximate  solution  over  a  finite  mesh 
of  elements  that  models  any  given  domain.  This  domain  can  include  virtually  any 
geometry  including  boundaries.  The  system  of  equations  is  assembled  from  all  dis¬ 
crete  finite  elements  of  the  domain,  and  solved  to  produce  a  solution  over  the  entire 
domain.  From  the  brief  description  given,  it  can  be  seen  that  the  fundamental  com¬ 
ponent  of  a  FE  solution  is  the  production  of  a  set  of  algebraic  equations  for  each 
element.  The  functional  form  of  the  coefficients  A  and  forcing  term  /  in  these  equa- 
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SAFE-T  Horizontal  vs.  Exact  Results 
For  Unit  Step  Function 


Figure  17.  Solid  Adaptive  Finite  Element  -  Transient  (SAFE-T)  is  used  to 
model  horizontal  displacement  due  to  an  instantaneous  Heaviside  change 
of  state.  Again,  roughly  92  percent  accuracy  is  achieved  with  2,764,800 
degrees  of  freedom.  At  time  steps  are  taken  at  0.0005  seconds.  This 
particular  calculation  was  done  on  a  Compaq  laptop  computer  with  an 
AMD-64  Athlon  processor  and  2  gigabytes  of  memory.  It  took  17.78  hours 
of  wall  time  with  3.05  hours  of  that  being  actual  CPU  time. 

tions  is  identical  for  each  element;  only  the  numerical  values  of  A  and  /  differ  from 
element  to  element.  Thus,  developing  a  FE  formulation  consists  of  developing  the 
functional  expression  for  A  and  /  for  a  master  set  of  element  equations,  which  can 
then,  like  a  master  template,  be  numerically  evaluated  over  and  over  again  for  each 
element  in  a  mesh  in  order  to  generate  the  assembled  system  of  algebraic  equations. 
The  theoretical  development  of  the  coefficients  can  be  found  in  great  detail  in  [Ref. 
33], 


a.  Commercial  FE  Development  Software 


The  principal  software  package,  Prophlex,  is  a  suite  of  developmental 
tool  for  creating  customized  finite  element  applications.  Solid  Adaptive  Finite  El¬ 
ement  -  Transient  (SAFE-T)  is  a  time  dependent  three  dimensional  finite  element 
tool  developed  by  the  author  using  Prophlex  for  analysis  of  wave  propagation  in  a 
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Rayleigh  wave 


Plane  S-wave 


Free  Surface 


P-wave 


Figure  18.  This  graphical  depiction  is  the  actual  element  by  element  de¬ 
formation  caused  by  a  Rayleigh  wave.  It  is  modeled  using  SAFE-T’s  nu¬ 
merical  output  and  postprocessed  using  the  Altair  Inc.  Hyperview  Post¬ 
processing  Suite.  As  the  wave  passes  an  observed  point,  all  three  waves 
(compressional,  shear,  and  Rayleigh)  can  be  seen.  The  values  of  A  and  fi 
for  this  particular  medium  is  equal  to  each  other  which  makes  v  0.25. 


solid  media.  Prophlex  is  applicable  to  any  system  or  process  that  can  be  mathemati¬ 
cally  formulated  into  a  system  of  second-order  PDEs.  The  FE  mesh  and  other  input 
data  for  the  SAFE-T  models  are  generated  using  Altair  Engineering  Inc.  developed 
Hypermesh  7.0. 

b.  Examples 


Figures  16  and  17  are  examples  of  SAFE-T’s  numerical  efforts  in  mod¬ 
eling  such  phenomena.  Notice  the  early  arrival  of  the  p-wave  in  the  SAFE-T  results  in 
figure  16  corresponds  to  the  exact  arrival  time  computed  analytically  with  an  error  of 
about  8  percent.  Figure  18  is  a  graphical  depiction  of  a  ’’scaled”  element  by  element 
deformation  of  a  Rayleigh  wave.  The  actual  amplitudes  are  on  the  order  of  microns. 
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2. 


Boundary  and  Initial  Conditions 


For  any  physical  problem  modeled  by  a  PDE,  many  solutions  are  possible 
given  the  variety  of  geometries  and  diverse  loads.  To  single  out  one  particular  solu¬ 
tion  requires  formulation  of  physically  realistic  boundary  and  initial  conditions  which 
together  make  the  problem  both  well-posed  and  useful  [Ref.  34],  The  FE  method 
requires  the  same  types  of  constraints  as  it,  too,  models  physical  problems  that  are 
themselves  well-posed.  As  such  three  types  of  boundary  conditions  are  needed  in 
order  to  model  the  clastic  half-space  for  this  problem. 

a.  Displacement  Ui  is  specified 

Ui  =  Ui  (11.105) 


b.  Normal  and  Tangential  Derivatives  of  displacements 
are  specified  ( applied  or  free  stress ) 

Ti:jh  =  ti(xi,t)  (11.106) 


c.  Normal  Derivative  of  displacement  is  zero 

5?  =  0  (11.107) 

on 

The  initial  conditions  for  this  problem  simply  state  that  the  medium  is  at  rest. 

3.  Time  Marching  Scheme 

The  time  marching  scheme  employed  is  based  on  the  Generalized  Newmark 
algorithm  (GNpj)  [Ref.  7].  This  time  marching  scheme  converges  to  the  analytic 
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Figure  19.  The  expansion  of  an  unknown  vector,  say  a,  will  be  taken  as 
a  second  degree  polynomial  with  known  values  for  an,  an,  and  an  at  the 
beginning  of  each  time  step  At. 

solution  as  the  time  step  shrinks.  Figure  19  is  a  graphical  depiction  of  how  succeeding 
solutions  are  calculated.  The  expansion  of  an  unknown  vector,  say  a,  will  be  taken  as 
a  second  degree  polynomial  with  known  values  for  an,  an,  and  an  at  the  beginning  of 
time  step  At.  The  Lax-Equivalence  Theorem  states  that  the  necessary  and  sufficient 
condition  for  a  numerical  method  to  be  convergent  are  consistency  and  stability  [Ref. 
35] .  Consistency  simple  means  that  as  time  intervals  between  calculating  solutions  is 
decreased,  the  truncation  error  between  the  numerical  solution  and  the  exact  solution 
approaches  zero.  Convergence  of  a  numerical  method  to  an  analytic  solution  implies 
that  the  numerical  method  is  consistent,  but  the  converse  is  not  true.  Consistency 
is  not  enough,  but  consistency  with  stability  is  enough.  Zero-stability  is  concerned 
with  the  stability  of  the  system  in  the  limit  as  the  time  intervals  between  calculating 
the  solution  shrinks  to  zero.  A  built-in  instability  exists  for  initial  value  problems 
even  in  the  limit  as  the  time  intervals  approach  zero  in  duration,  but  is  mitigated 
by  ensuring  that  the  roots  of  the  characteristic  equation,  or  root  condition,  of  the 
numerical  method  have  absolute  magnitudes  that  are  less  than  or  equal  to  unity. 
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According  to  Isaacson  and  Keller  [Ref.  36],  fact  that  the  highest  order  of  accuracy 
that  can  be  expected  from  a  k- step  method  is  2k,  so  a  single  step,  k— 1,  numerical 
method  can  be  at  most  second  order  accurate.  Use  of  a  single  step  method  conserves 
computer  resources  by  calculating  solutions  with  one  pass  of  the  computer  processor. 
While  it  may  be  true  that  in  seeking  higher  order  the  consistency  condition  is  well 
satisfied,  attempting  to  satisfy  the  condition  for  zero-stability  becomes  impossible. 
This  barrier  was  first  investigated  by  Germund  Dahlquist  [Ref.  37]  and  expounded 
upon  by  Lambert  [Ref.  38].  Because  of  this  barrier,  the  second-order  single-step 
variant  of  GNpj  called  the  Newmark  Beta  Method  gives  the  highest  order  achievable 
in  a  single  step.  It  is  derived  by  making  approximations  for  velocity  and  acceleration 
terms  that  come  from  third  order  Taylor  Series  expansion  at  the  tn+l  time  step. 


u 


7+1  =  u?  +  At 


du'l 

dt 


At2  d2vA  At3  d3u™ 

+  ^^t  +  iri^  +  0(h) 


+  ... 


(11.108) 


The  third  derivative  term  is  approximated  by  a  backwards  difference  scheme  leading 
to 

3vn  At2  d2vn  ( d2vn~^3  d2vn\ 

u?+1  =  u?  +  A f wrxf  +  (3A 12  - )  +  0(At4)  +  . . .  (11.109) 


dt 


2  dt 2 


dt 2 


u 


dvn  /I  \  d2vn  d2vn+1 

n+1  =u?  +  A +  At2  ( -  -  (3 )  ^  +  I3A t2^r-  +  0(At4)  +  . . .  (11.110) 


dt 


dt2 


dt 2 


g‘2un+1 

and  now  solving  for  the  at 2  term  gives 


d2u™+l 
dt 2 


(3  At 


1 


(3 At  dt 


2  (3 


d2u " 
~d¥ 


+  0(At2)  +  . . .  (11.111) 


and  similarly  for  the  first  derivative  term  we  have 


dK+1 

dt 


7 


(3At 


(<+1  -  <) 


7 


[u?+1  -u?)  -  [-  -1 

\d  / 


dun{ 

~dt 


-  At 


1_ 
2  (3 


-  1 


d2u 
dt 2 


+  0(At2)  + 


(11.112) 

The  correct  representations  for  derivative  terms  are  now  used  to  perform  time  march¬ 
ing  in  a  one-step  scheme.  Thus,  we  accomplish  a  more  economic  use  of  computer 
resources  as  each  time  step  can  be  solved  in  one  pass  of  the  computer  processor. 
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Surface  Response  1m  away 
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Figure  20.  SAFE-T  vertical  displacement  results  compared  to  the  arbitrary 
source  solution  derived  by  discrete  integral  convolution.  In  this  model  an 
extended  mesh  is  used  to  prevent  unwanted  body  waves  from  interfering 
with  the  result. 

A  transient  numerical  approach  has  been  developed  that  incorporates  initial 
and  boundary  conditions  that  makes  the  mathematical  model  both  realistic  and  use¬ 
ful.  By  discrete  integral  convolution,  an  analytic  benchmark  has  been  computed  to 
verify  the  accuracy  of  the  FE  method.  Figure  20  is  a  comparison  of  SAFE-T’s  vertical 
displacement  results  to  the  analytic  benchmark.  This  result  is  accomplished  through 
the  time- marching  scheme  we  have  just  discussed.  However,  this  is  only  half  of  the 
story.  The  other  half  centers  about  the  difficult  task  of  truncating  the  computational 
domain  in  such  a  way  that  it  models  an  infinite  half-space.  We  accomplish  this  by 
introducing  a  truncated  elastic  media  via  perfectly  matched  layers. 
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III. 


TRUNCATED  ELASTIC  MEDIA 


A.  INTRODUCTION 

A  new  unsplit  3D  time  dependent  elastic  PML  PDE  will  be  derived  and  con¬ 
verted  to  its  weak  (Galerkin)  form  for  implementation  into  the  finite  element  method. 
We  have  so  far  offered  the  finite  element  method  as  a  numerical  approach  to  the  ap¬ 
proximate  solution  of  the  elastic  PDE  with  appropriate  boundary  conditions.  This 
is  well  established  as  a  reliable  choice  for  problems  that  are  finite.  Infinite  or  semi¬ 
infinite  problems  are  impossible  without  some  way  to  absorb  undesired  reflections. 
Figure  21  is  an  example  of  the  devastating  effects  of  body  wave  reflections  that  makes 
reliable  wave  propagation  modeling  impossible.  The  graph  represents  the  displace¬ 
ment  history  of  a  point  on  the  free  surface  of  the  elastic  half  space.  The  non-PML 
response  matches  the  PML  response  up  to  40  At  time  steps,  but  soon  reflection  of 
body  and  surface  waves  corrupt  the  time  histories.  The  reflections  distorts  the  am- 


Figure  21.  The  devastating  effects  of  body  waves  make  reliable  wave  prop¬ 
agation  modeling  impossible. 


43 


plitude  of  the  surface  wave  and  renders  it  useless  as  a  scatterer.  The  purpose  of  this 
chapter  is  to  show  that  all  reflections  of  body  waves  as  well  as  reflected  surface  waves 
can  be  suppressed  and  absorbed  by  a  transient  PML  boundary. 

The  objective  of  any  PML  method  is  to  construct  a  new  wave  equation  that 
causes  waves  to  decay  exponentially  as  they  traverse  the  PML  [Ref.  39].  This  is 
accomplished  most  effectively  by  a  complex  change  of  variable  for  xi  as 

xi — +  MQdt)  (nL1) 

where  fi(x)  >  0  is  the  absorption  function  and  a  is  the  location  of  the  PML  interface. 
By  replacing  X/  in 

um(xl,u)  =  Amei^k^  (III.2) 

we  have 

Um(xh  u)  =  (III.3) 

This  transformation  has  the  desired  effect  of  introducing  a  purely  real  exponential 
term  into  the  expression  which  acts  to  decay  wave  amplitude  while  the  eigenvalues  of 
the  contained  (computational  domain)  remain  unchanged  [Ref.  40].  If  we  let  a(xi) 
equal  the  non-zero  real  quantity  used  to  suppress  the  wave  such  that 

fci  rxi 

a(Xl)  =  -*  /  /,(£)#  (HL4) 

LO  Ja 

then  the  damping  will  be  in  the  direction  of  ay.  For  example,  suppose  one  wanted 
to  decrease  a  plane  wave  traveling  in  the  aq-direction.  A  new  governing  equation  is 
required  whose  plane  wave  solution  would  have  the  form 

Umix^u)  =  Ame^xikl~^-a^  (III. 5) 

where  a(xi)  is  a  value  that  is  nonzero  only  in  the  PML  region  and  can  be  adjusted 
to  produce  a  desired  decay  rate.  The  purpose  of  this  complex  coordinate  stretching 
variable  is  to  alter  a  wave’s  behavior  as  it  traverses  the  PML  region.  Another  way 
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to  derive  the  complex  coordinate  stretching  variable  was  introduced  by  Chew  and 
Weedon  [Ref.  16]  as 

•0  f  /•}(«•  (j  =  1,2,3)  (III.6) 

Jo 

Fj{x)  =  l-ifj(x)  (III.7) 

Here  /  is  an  attenuation  factor  in  the  PML  region  and  is  zero  within  the  computa¬ 
tionally  significant  domain.  Thus, 


j)  =  I J  /-.'/(sX- 

Jo 


(III.8) 


is  used  to  produce 


Xj  =  Xj  —  a{Xj/ 


(111.9) 


where  the  stretching  coefficient  along  the  prescribed  axis  Xj  is  a  complex  number. 
The  resulting  independent  variable  causes  damping  of  wave  fronts  propagating  in 
a  prescribed  direction,  and  is  very  efficient  for  wave  absorption  at  the  boundary 
of  a  numerical  model  [Ref.  41].  The  complex  coordinate  stretching  function,  F, 
is  continuous,  and  therefore,  the  stretch  coordinate  is  smooth  and  the  Fundamental 
Theorem  of  Calculus  can  be  applied.  It  can  be  shown  that  if  we  differentiate  equation 
III. 6  with  respect  to  Xj, 

f)'V 

(III.  10) 


v3 


which  implies  that 


_d_ 

dxj 


]__d_ 

FjSij 


(in.11) 


We  now  have  a  new  differential  operator  for  the  governing  equations  that  incorporate 
the  PML  properties  of  the  media  through  its  differential  operator.  The  PDE  becomes 


1  <9r, 


Fj  9xj 


1-3  2- 

=  ~pU  U  i, 


(III. 12) 


and  the  constitutive  equations  are 


7~ij  Cijkl^kl  ? 


(III. 13) 
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f  j  (x) 


Figure  22.  Linear  and  nonlinear  damping  functions  used  in  the  PML. 


and 


6ij  ~  2 


1  1 


(III.  14) 


where  F  has  a  complex  value  only  in  the  absorbing  PML  region  and  is  otherwise 
unity. 


At  the  interface  of  the  computational  domain  and  the  PML,  the  wave  equations 
are  identical  so  that  any  propagating  wave  will  pass  through  the  interface  without 
generating  reflected  waves.  In  other  words,  there  is  an  impedance  match  at  the 
boundary  of  the  unsealed  computational  domain  and  the  PML.  Care  must  be  taken 
when  choosing  an  appropriate  Fj.  It  must  be  complex  as  defined  in  equation  III.  7 
with  an  imaginary  part  related  to  the  desired  wave  attenuation  [Ref.  42].  It  is  also 
desirable  to  have  the  imaginary  part  of  equation  III. 7  increase  gradually  relative  to  its 
position  in  the  PML.  This  will  provide  zero  attenuation  at  the  interface  yet  gradually 
increase  the  attenuation  as  the  wave  travels  in  the  direction  of  the  PML  x3  coordinate. 
Choosing  f(x)  to  be  linear  or  nonlinear  in  the  PML  region  are  possible  candidates 
for  fulfilling  the  attenuation  requirement.  Figure  22  is  a  graph  of  an  example  of  both 
the  linear  and  nonlinear  functions  fj(x).  In  the  linear  case,  0  <  f3  <  1 


fj(x)  = 


t xj  —  P  M  LStart 
,  LpML 


if  the  plane  wave  is  in  the  PML 


if  the  plane  wave  is  in  the  computational  area 


(III. 15) 


In  the  nonlinear  case,  0  <  f3  <  1 
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Comparison  of  Stretching  Functions 


Figure  23.  Comparison  of  linear  7  =  1,  quadratic  7  =  2,  and  cubic  7  =  3 
damping  functions  used  within  the  PML  region  of  the  FE  model.  All 
three  functions  whether  linear  or  non-linear  provide  excellent  absorbtion 
with  almost  identical  results. 

(  ( Lj-PMLstarty  if  the  lane  wave  ig  in  the  PML 

fj{x)  =  l  V  Lpml  J  (III.  16) 

1  0  if  the  plane  wave  is  in  the  computational  area 

where  PALLStart  is  where  the  PML  region  begins,  Lpml  is  the  total  length  of  the 
PML  region,  7  is  a  non-linearity  constant,  a0  is  a  damping  constant,  and  lxj  is  the 
length  from  the  interface  to  the  plane  wave  front.  Figure  23  is  a  comparison  of  the 
effects  of  using  linear,  quadratic,  and  cubic  functions  as  the  damping  functions  within 
the  PML  region.  All  three  functions  whether  linear  or  non-linear  provide  excellent 
absorbtion  with  almost  identical  results.  Figure  24  analyzes  the  sensitivity  of  the  pml 
to  the  damping  constant,  «o- 
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PML  Alpha  Analysis 


Figure  24.  Analysis  of  the  PML  to  determine  the  sensitivity  to  the  damping 
constant,  a0.  The  location  of  measurement  is  along  the  computational 
domain/PML  interface.  The  bump  circled  is  highlighting  the  reflected 
amplitude  of  the  incident  wave  upon  the  boundary  because  of  a  stiffening 
of  the  interface,  a o  set  to  2  provides  the  maximum  absorption  with  minimal 
interface  stiffness  and  incident  reflection. 

B.  3D  TRANSIENT  PML  EQUATION  DERIVATIONS 


Although  PMLs  are  more  complex  to  implement,  especially  in  three  dimen¬ 
sions,  the  methodology  for  formulating  their  equations  are  roughly  the  same.  PMLs 
require  only  a  finite  number  of  nodes  per  wavelength.  Between  4  and  8  nodes  per 
wave  length  provide  an  excellent  absorption  of  body  waves,  and  do  not  lose  efficiency 
at  shallow  angles.  Most  notably,  PMLs  have  been  shown  to  be  very  effective  with 
surface  waves.  The  reflections  at  the  boundaries  can  be  made  arbitrarily  small  by 
increasing  the  thickness  of  the  PML  layer  at  the  cost  of  additional  computation  [Ref. 
41].  This  cost,  however,  is  usually  well  worth  the  extra  resources  required.  With  that 
in  mind,  we  adopt  a  perfectly  matched  medium  undergoing  time-harmonic  motion  in 
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the  absence  of  body  forces,  whose  governing  equations  III.  12,  III.  13,  and  III.  14,  are 
defined  by  the  following: 


Fl(Xl)  F2(X2)  F3(x3)  -  2  cW  x£w  xj=w  X- 

- p  7  x - rij,j  =  -U  pF1{Xi)F2{X2)F3{x3)Ui 

rj{Xj) 

Tij  Cijkl^kl 

1 


eij 


(111. 17) 

(111. 18) 

(111. 19) 


1  1 

2  [fmu,j + 

The  stretching  function,  F^Xi),  must  possess  special  properties.  It  must  be  unity  in 
the  computational  domain  and  complex  otherwise  with  a  real  component  that  damp¬ 
ens  evanescent  waves  and  a  complex  component  that  dampens  propagating  waves 
in  the  PML.  A  detailed  analysis  of  stretching  functions  can  be  found  in  [Ref.  11], 
Choosing  the  stretching  function  to  be  of  the  form 

Csfi{Xi)' 


Fi{xi )  —  [1  +  fi(%i )]  + 


IU) 


(III. 20) 


where  fe  and  fp  are  evanescent  and  propagating  damping  functions  respectively,  and 
cs  is  the  shear  wave  phase  speed-  used  as  the  reference  speed  in  the  clastic  medium, 
yields  stretch  matrices  A,  A  and  A.  These  three  matrices  contain  all  the  coordinate 
stretching  information  of  the  perfectly  matched  medium.  When  in  the  computational 
domain  they  are  zero,  (A  and  A),  or  identity  (A),  but  in  the  perfectly  matched 
medium  they  attenuate  both  evanescent  and  propagating  waves. 


A  = 


( 


\ 


(l  +  /2e)(l  +  /|)  0 

0  (i  +  /i)(i  +  /r; 

0  0 


0 

0 
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(i  +  /f)(i  +  /f)  j 


A  = 


cs[fi  (1  +  /I)  +  /I  (1  +  /!)]  0 

0  cs[/f(l  +  /f)  +  /f(l  +  /|)] 
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0 


0 
0 

cs[fl(  1  +  /I)  +  f$0-  +  ft)] 


A  = 


(  r2  fP  fP 
lsJ2 J3 
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c2JU 'f  0 
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The  scaled  or  stretched  governing  partial  differential  equation  in  the  frequency  domain 
becomes 


+ 


A jmTij,m 

iuj 


A-jrriTij.m, 


UJ “ 


—  — c o2p  [1  +  /f]  + 


Csfl 

iuj 


[1  +  /I]  + 


iuj 


[i + m  + 


c_Jt 

iuj 


Uj 


(III. 21) 


and  contains  within  it,  all  the  structure  needed  to  implement  a  PML  along  any  of 
the  three  coordinate  directions  in  a  rectangular  system.  Because  the  shear  speed  is 
chosen  as  the  reference  speed,  shear  and  surface  waves  have  near  perfect  absorbtion. 
Compressional  waves  which  travel  at  greater  phase  speeds  have  minor  reflections  that 
can  be  made  arbitrarily  small  by  adjusting  the  thickness  of  the  PML. 


1.  Integral  and  Inversion  Techniques 

a.  Frequency  Inversion 


PML  equations  for  the  frequency  domain  are  given  in  equation  III. 21. 
In  order  to  use  the  PML  method  in  the  time  domain,  a  transformation  must  occur. 
Since  multiplication  or  division  by  the  factor  iuo  in  the  frequency  domain  is  equiva¬ 
lent  to  differentiation  or  integration  respectively  in  the  time  domain,  the  equations 
are  transformed  into  their  time-dependent  counterparts  by  application  of  equation 
11.36,  the  Fourier  inversion  formula  [Ref.  14],  The  application  of  the  inverse  trans¬ 
form  assumes  that  r  is  zero  when  ui  =  0.  With  this  assumption,  equation  III. 21  is 
transformed  from  the  frequency  domain  PML  equation  of  motion  to  the  time  domain 
PML  equation  of  motion 


A jmTijjm  A-jm  /  7 A-jr 


t  rC 


o  J  o 


7~ij,md^dC 


=  Miii  +  Dili  +  Kut  +  L  I  utd^ 

Jo 


(III.22) 


where  following  the  convention  used  in  [Ref.  14]  in  2  dimensions  and  applying  it 
to  3  dimensions  yields  a  mass  term,  M,  a  damping  term,  H,  a  stiffness  term,  K, 
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and  a  time-integral  term,  L.  These  are  somewhat  unconventional  from  a  continuum 
mechanics  view,  but  naturally  arise  in  a  time-domain  implementation  for  a  PML, 
according  to  Zhao  (1996),  when  field-splitting  is  avoided  [Ref.  43]. 

M  =  p(l  +  /f)(l  +  /|)(1  +  /I) 

n  =  pc.[/f(i  +  /2e)(i + /i)  +  m  +  /r)(i  +  m  +  /f(i + /r)(i  +  m 

K  =  pc2s[f2f f  (1  +  ft)  +  /f/f  (1  +  /I)  +  /i/f(l  +  /!)] 

£  =  pclfUUl 

By  making  a  substitution  for  the  integrals,  a  3-D  transient  clastic  wave  equation 
containing  the  PML  parameters  can  be  written  in  the  complete  form 

in  T~ij.ni  "P  A_ym  ^  ij. in  HUi  T  Dili  DU{  T  DU i  (III. 23) 

where 

T/j.m  j  j •  ^ij,m  J  j  ■  D /  J Uidf  (III. 24) 

b.  Weak  Formulation 


Construction  of  a  weak  formulation  of  III. 23  begins  with  factoring  di¬ 
vergence  terms  and  consolidating  them  on  the  LHS.  This  leads  to 

(A  jmJij  +  A  jm&ij  +  A  jm^ij),m  =  MHz  +  Dili  +  KUi  +  LTj  (III. 25) 


Proceeding  in  the  usual  manner,  we  multiply  III. 25  by  a  test  function  Vi  and  integrate 
over  the  domain,  Q: 

[(A T  A jm^ij  +  A jmi&ij),m]vidQ  =  Ja(Mui  +  Dili  +  Kui  +  LTfjVidD  (III.26) 

By  Green’s  Identity  (integrating  by  parts)  we  expand  the  terms  on  the  LHS,  and  are 
left  with 


j  (A jmTij  A  jm^ij  ~\~  A  jm^ 

=  [  ( Mui  +  Dili  +  Km  +  LTAvidn  (III.27) 

Jn 
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Finally,  applying  the  divergence  theorem,  we  derive  a  weak  formulation  of  the  3-D 
transient  elastic  wave  partial  differential  equation  that  incorporates  the  PML  param¬ 
eters.  Thus, 


/  T  T  )uj]'hm(ir  /  (AjmTij  T  T  Aj , n T jj  ) Vi  r n 0 1 1 

Jr  Jn 

=  [  ( Mut  +  Dili  +  Kiii  +  LTi)vidQ  (III. 28) 

Jn 

where  the  first  term  in  equation  III. 28  is  a  boundary  surface  integral. 


2.  Galerkin  Surface  Integral 


Consider  the  surface  integral  equation  III. 28. 

/  [(AjjtjTjj  T  A jm&ij  -(-  Ajm^f ij')Vi]nmdT'  (III. 29) 

j  r 

The  integration  is  only  over  boundary  surface  elements.  For  3-D  brick  elements, 
integration  is  over  at  most  three  faces  of  any  brick.  Because  of  the  presence  of  the 
PML,  this  really  leads  to  three  possible  conditions  that  may  be  prescribed  on  an 
element  face  as  seen  in  Figure  25.  Equation  III. 28  may  be  expanded  to  include  the 
possibilities: 


“1“  Ajm^ij  ^ D  “I-  ^ N 

(III. 30) 

where 

To  =  J 

(III. 31) 

Tn  =  J 

+  AjfyJ&ij  +  Ajm  Tlm  C?r , 

Fn 

(III. 32) 

and 

rT  1  “1“  A-jm^ij  “1“  Ajm^ij 

JYt 

represent  Dirichlet,  Neumann,  and  Stress  conditions,  respectively. 

(III. 33) 
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Stress  Boundary 


Homogeneous  Dirichlet 


Figure  25.  Iso  view  of  a  FE  model  of  a  solid  with  Stress  and  Dirichlet 
boundaries.  Any  one  of  two  conditions  may  be  present  on  an  element 
face. 


a.  Condition:  Dirichlet 


The  Dirichlet  boundary  condition  assigns  specific  displacement  (in)  val¬ 
ues  on  the  boundary.  Thus, 


Ui  t/j. 


(III. 34) 


However,  this  representation  cannot  be  directly  applied  to  the  Dirichlet  integral  on 
the  RHS  of  the  expanded  boundary  integral  (III. 30).  This  is  because  the  boundary 
integral  is  in  terms  of  normal  derivatives  of  stress  instead  of  tp.  This  condition  can 
be  circumvented  by  developing  a  penalty  parameter  formulation  of  the  integral.  It  is 
accomplished  by  specifying  that  the  difference  between  Ui  and  U,  be  a  small  e  quantity 
such  that 

Ui  Ui  c(A jmTij  T  A jm^ij  T  Ajm^ij)nm  (III. 35) 

-(Ui  Ui)  (A jmTij  T  A jm^ij  T  A jrn^ij)^'m  (III. 36) 
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Figure  26.  Dual  iso  views  of  the  FE  model  of  a  solid  homogeneous  cube. 
The  boundaries  of  the  cube  are  Dirichlet  with  a  free  surface  directed  up¬ 
ward.  Notice  the  symmetry  of  wave  propagation.  The  source  is  a  Gaussian 
ignited  at  the  center  of  the  domain.  Wave  propagation  is  governed  by  the 
medium. 

where  e  is  the  penalty  parameter  [Ref.  7].  Equation  III. 36  can  be  substituted  directly 
into  III. 31  to  give  us  a  new  representation  for  T /> 

l  [(A  jmXij  A-jm  tf  j-j )  U7]  h;ndh  J  UjVj(]Y  j  UjVjflT'  ( 1 1 1 . 3  7  [) 

The  value  of  e  is  chosen  to  be  extremely  small.  This  is  to  make  sure  that  when  III. 37 
is  applied  to  the  final  assembly  matrix  the  two  integrals  on  the  RHS  of  III. 37  will 
dominate  all  other  algebraic  terms  in  the  equation.  The  e  domination  will  leave,  in 
effect,  only  the  two  4  integrals  in  which  case  the  coefficients  cancel  each  other  out. 
This  leaves  us  with  the  condition  defined  by  equation  III. 34.  Typically,  the  guidelines 
for  choosing  e  is  somewhere  around  eight  to  ten  significant  digits  [Ref.  44], 
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b. 


Condition:  Neumann 


The  Neumann  boundary  condition  might  be  used  along  the  exterior 
boundaries  of  the  PML  but  if  universally  applied  could  lead  to  a  non-unique  solution 
of  the  governing  equation.  Some  portion  of  the  PML  or  free  boundaries  must  be 
assigned  Dirichlet  conditions  to  ensure  uniqueness  of  the  numerical  results. 

c.  Condition:  Stress 


The  stress  condition  applies  a  specific  value  to  normal  derivatives  of 
displacements  at  the  bounded  surface  which  may  constitute  a  free  surface  when  equal 
to  zero  or  an  applied  known  stress  when  inhomogeneous.  We  derive  it  by  substituting 


( A-jmTij  T  A jm^ij  T  A-jm^ij )^m 
into  the  RHS  of  111.33  yielding 


ti(xi,t),  applied  stress 
0  stress  free 


(111.38) 


(111.39) 


where  U(xi,t )  specifies  a  specific  time  dependent  source  of  stress. 


3.  Boundary  Conditions 

The  Galerkin  surface  integral  of  this  model  defines  two  surface  conditions: 
stress,  and  displacement.  Because  the  free  stress  condition  contributes  nothing  to  the 
RHS  of  111.30,  stress  (where  a  non-zero  forcing  function  is  applied)  and  displacement 
are  the  only  two  boundary  conditions  that  need  be  specifically  calculated.  The  ge¬ 
ometry  of  the  half-space  is  very  simple  and  requires  no  special  treatment.  As  such, 
a  fixed  Dirichlet  condition  is  suitable  for  the  outer  boundary  of  the  PML  boot  which 
completely  surrounds  the  computational  domain  except  for  the  free  surface  of  the 
half-space  (see  figure  25).  The  free  surface  source  or  stress  traction  requires  more 
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Figure  27.  Top  view  of  a  FE  model  of  a  solid  cube  with  dirichlet  boundaries. 
Notice  the  symmetry  of  wave  propagation.  The  surface  wave  will  impinge 
upon  the  boundary  and  return  180  degrees  out  of  phase  and  travel  back 
across  the  domain.  This  totally  unphysical  mathematical  phenomenon  is 
exactly  what  the  PML  is  designed  to  solve.  We  say  unphysical  only  in  the 
sense  that  in  an  infinite  half-space,  waves  travel  outward  and  never  return 
to  the  originating  source. 
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Figure  28.  Top  view  of  a  FE  model  of  a  solid  cube  with  a  PML  boot 
truncated  by  Dirichlet  boundaries.  Again,  notice  the  symmetry  of  wave 
propagation.  The  surface  wave  will  impinge  upon  the  boundary  and  be 
totally  absorbed.  This  is  exactly  what  the  PML  is  designed  to  do.  Waves 
travel  outward  and  never  return  to  the  originating  source. 
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attention.  The  forcing  function  of  the  half-space  problem  is  a  point  load  directed 
downward.  Figure  27  shows  the  reaction  of  waves  excited  by  a  point  source  imping¬ 
ing  upon  a  Dirichlet  boundary  without  a  PML  present.  It  represents  one  shaker  that 
exerts  a  vertical  traction  only.  This  source  is  always  placed  within  the  computational 
domain  away  from  the  PML  so  that  the  stretch  tensors  are  always  zero  or  identity. 
Recalling  the  weak  form  of  the  elastic  equation  with  PML  parameters  added  we  see 
that  the  only  forcing  term  is  the  boundary  surface  integral,  Ft,  which  is  the  normal 
component  of  the  surface  stress  tensor.  Figure  28  shows  the  reaction  of  waves  excited 
by  a  point  source  impinging  upon  a  Dirichlet  boundary  with  a  PML  present.  Note 
that  to  eliminate  having  a  rectangular  wavefront  in  figures  27  and  28  which  is  the 
result  of  grid  dispersion,  a  finer  mesh  is  used  to  produce  a  more  cylindrical  shape. 
A  rectangular  shape  occurs  when  course  meshes  are  used  to  model  relatively  fast 
wave  phenomena.  The  graphs  above  clearly  demonstrates  the  efficacy  of  the  PML 
boundary. 


4.  Time  Integration 


The  single-step  recurrence  relations  derived  in  Chapter  II,  equations  11.110, 
11.111,  and  11.112,  are  applied  to  the  PML  PDE,  equation  III.28.  By  substituting  the 
approximations  for  the  first,  and  second  time  derivatives  into  the  weak  form  of  the 
PDE  and  collecting  unknown  terms  (primarily  those  at  time  tn+1)  on  the  LHS  and 
collecting  terms  at  time  tn  and  earlier  on  the  RHS  yields  a  discrete  approximation 
to  the  solution  of  all  primary  variables.  Note  that  stress  and  displacement  boundary 
conditions  are  both  known  at  tn+1  and,  therefore,  are  placed  on  the  RHS.  Thus, 
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(III. 40) 
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is  a  diagonal  stretching  matrix.  The  trapezoidal  rule  [Ref.  45]  is  used  to  approximate 
integrals  with  global  truncation  errors  of  0(Af2): 
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(111. 42) 

(111. 43) 

(111. 44) 


i=i 


C.  STRAIN  AND  THE  PML 

At  this  point  it  is  necessary  to  calculate  the  strain  of  a  system  in  the  PML. 
In  order  to  mathematically  capture  the  properties  of  the  perfectly  matched  damping 
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media,  a  new  stretching  matrix  will  be  defined.  This  matrix  will  allow  the  PML 
properties  of  the  media  to  be  combined  with  the  classic  strain  equation  in  order  to 
express  a  new  strain  equation.  We  begin  with  a  diagonal  matrix,  say  S,  which  consists 
of  the  reciprocals  of  the  stretching  functions  defined  in  III. 20 
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1  1/ Fi(xi)  0  0  ^ 

0  l/F2(x2)  0 

v  0  0  1/F3(x3)  f 


(III. 45) 


Matrix  S  will  be  the  mathematical  vehicle  used  to  transfer  perfectly  matched  layer 
properties  into  the  classic  strain  equation.  Furthermore,  the  summation  convention 
will  be  abandoned  in  certain  cases  below  and  in  some  cases  matrix  to  matrix  products 
will  be  found  by  multiplying  term  by  term  elements  in  each  matrix  to  produce  a 
product  matrix.  For  example,  the  i,  j  component  of  C  =  AB  will  become  AjjBjj  = 
C ij.  With  that  in  mind;  using  matrix  S  above,  we  can  rewrite  the  classic  strain 
equation  to  yield 

%  =  (III. 46) 

Standard  left  and  right  multiplication  by  the  inverse  of  the  stretching  matrix  S  yields, 


S-^S"1  =  ^(S^fqjSS-1  +  S_1S%iS-1)  (III. 47) 

S-^-S-1  =  ^(S-1^-  +  uj.  S’1).  (III. 48) 


Basu  and  Chopra(2003)  [Ref.  14]  make  ample  use  of  this  technique  for  ID  and 
2D  strain  manipulations.  Here,  for  the  3D  case,  we  make  use  of  the  technique  to 
isolate  the  strain  terms  by  the  factor  iuj  with  the  substitution  of  the  stretch  function 
F.fxt),  defined  in  equation  III. 20.  By  placing  this  into  equation  III. 48  and  abandoning 
summation  over  double  indices  in  the  expression  below,  we  have 
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Algebraic  manipulations  of  the  strain  and  displacement  terms  are  grouped  in  powers 
of  the  frequency  domain  variables  which  await  transformation  to  the  time  domain 
(again,  there  is  no  summation  over  the  indices).  Thus, 


[i +  /f][i+/;]  + 


c  (  fP\  1  _i_  fe 1  4.  f'P\ l  4.  f?])  2  rp  fp\ 
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(111.50) 


As  pointed  out  by  Basu  and  Chopra  [Ref.  11],  one  can  make  use  of  the  fact  that 
transformations  from  the  frequency  domain  to  the  time  domain  are  simple  when 
given  terms  that  are  multiples  of  iu.  Further,  if  stretch  matrices  B,  B  and  B  are 
defined 

'  (i  +  /f)(i  +  ft)  (i  +  /d (i  +  ft)  (i  +  /f ) (i  +  /i)  N 

B  =  (i  +  /!)(i  +  /f)  (1  +  /JK1  +  /I)  (1  +  /2K1  +  /3) 

y  (1  +  /|)(1  +  fl)  (l  +  /|)(l  +  /l)  (1  +  /|)(1  +  /|)  ) 
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As  remarked  above,  we  will  abandon  the  summation  convention  in  favor  of  an  element 
by  element  product  convention.  The  strain  is  discretized  in  time  and  approximated 
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using  the  trapezoidal  rule  for  the  integral  terms  such  that 
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(111. 53) 

(111. 54) 


fiui,j  +  /f%i]  At  J2  [fiui,j  +  f]u\i 

i=i 

Using  the  above  approximations  for  the  integral  terms  produces  an  equation  for  strain 
that  incorporates  the  PML  parameters.  The  fJn+X  time  level  strain  is 
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(III. 55) 


(III. 56) 


and 


Fi(xt)  =  ([1  +  /?(*,)]  +  ^/f fe))  .  (HI- 57) 

Equation  III. 55  emphasizes  a  change  in  material  parameters  when  used  to  derive 
stress.  Since  stress  is  proportional  to  strain,  the  perfectly  matched  media  can  thus  be 
interpreted  as  a  medium  which  exhibits  inhomogeneous  elastic  properties.  Note  that 
outside  the  PML  the  stretching  functions  ff{xj)  and  /f  (ay)  are  identically  zero  and 
the  above  strain  tensor  collapses  to  its  classic  form.  When  the  above  strain  is  placed 
into  the  stress-strain  equation,  we  have 
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(III. 58) 


where  the  additive  term  x%j  is  defined  as 
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Hankel  Plot  of  a  Surface  Particle  1m  Away 


Figure  29.  This  Hankel  plot  is  of  a  particle  on  the  surface  of  a  half-space 
r  distance  away  from  the  source.  It  clearly  shows  that  as  the  wave  moves 
across  the  surface,  a  particle  moves  in  a  circular  or  elliptic  pattern  in  the 
direction  of  propagation  which  is  characteristic  of  Rayleigh  waves. 

The  current  stress,  t/]+1,  can  now  be  used  in  the  weak  Galerkin  formulation.  Notice 
again  that  inside  the  computational  domain,  the  value  of  Xij  is  zero  as  expected  since 
inside  the  computational  domain  there  is  no  artificial  damping.  If  the  value  of  Xij  were 
not  zero  outside  the  PML  region  of  the  model,  damping  or  exponential  growth  would 
be  present  that  might  cause  unrealistic  growth  or  decay  in  regions  of  the  domain  of 
interest. 
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D.  FINAL  GALERKIN  FORM 


The  Galerkin  (weak)  form  of  the  problem  is  arranged  by  placing  all  terms 
with  implicit  components  on  the  LHS  and  explicit  components  on  the  RHS.  This 
form  of  the  equation  provides  all  the  information  needed  to  set  up  the  FE  model  in 
SAFE-T  using  the  Prophlcx  kernel.  Figures  29  and  2  are  results  obtained  by  SAFE-T 
using  ideal  material  properties.  The  problem  is  idealized  because  in  this  instance  the 
material  properties  A  and  //  are  set  equal  to  each  other. 
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hi\  =  tc\M  T  tc^D  (III. 61) 

k2  =  tc2M  +  tc$D  (III. 62) 

^3  =  tc3M  +  tc$D  (III. 63) 


Simplifying  and  collecting  terms  yields  the  final  compact  form  of  the  equation: 
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Surface  Response  1m  away 


Figure  30.  A  demonstration  of  the  remarkable  effects  of  the  PML  shows 
that,  although  not  identical  to  the  analytic  solution,  the  result  is  very 
close  to  the  extended  mesh  solution.  The  graph  shows  that  the  PML  has 
a  slight  effect  on  the  results  obtained  within  the  computational  domain. 
This  is  due,  in  part,  to  the  fact  that  the  P-wave  is  not  perfectly  matched. 


The  final  Galerkin  form  is  implemented  into  SAFE-T  by  way  of  Prophlex  through 
FORTRAN  and  C++  subroutines. 


E.  SAFE-T  RESULTS:  PML  EFFECTIVENESS 


The  true  benefit  of  the  PML  method  is  its  ability  to  allow  the  domain  to  be 
smaller  than  using  an  extended  mesh.  Extended  meshes  are  costly.  For  example, 
if  one  wanted  to  use  an  extended  mesh  to  model  elastic  wave  phenomena  in  sand, 
you  would  need  a  domain  greater  than  800  meters  to  model  1  second  of  Rayleigh 
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activity.  For  ID  problems  that  would  require  6400  elements.  In  3D,  it  is  triple  that 
amount.  The  reason  is  that  the  P-wave  travels  about  1600  meters  per  second  in  sand 
while  the  Rayleigh  wave  moves  only  95  meters  per  second.  The  remaining  hgures  and 
tables  of  this  chapter  are  an  analysis  of  the  thickness  of  a  PML  truncating  an  elastic 
computational  domain. 


x10  3  PML  Thickness  Analysis 


Time  In  Seconds 


Figure  31.  SAFE-T  PML  Thickness  Analysis  (Vertical  Displacements) 
shows  that  a  PML  that  is  3m  in  depth  is  very  similar  to  a  PML  5m  in 
depth.  This  resource  savings  is  not  trivial.  The  source  is  the  derivative  of 
a  Gaussian.  It  has  a  dominant  frequency  of  599.675  Hz.  The  characteristic 
wavelength  is  1  meter  with  a  mesh  density  of  8  nodes  per  meter  before 
refinement.  The  At  time  step  of  0.0005  is  well  within  the  CFL  stability 
criteria. 
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Figure  32.  SAFE-T  PML  Thickness  Analysis  of  (Total  Nodal  Displacement 
Magnitudes).  The  bump  at  about  .04  seconds  shows  the  early  arrival  of 
reflected  waves  from  the  fixed  boundary  truncating  the  PML.  Again,  a 
PML  that  is  3m  in  depth  is  very  similar  to  a  PML  5m  in  depth.  The 
characteristic  wavelength  is  1  meter  with  a  mesh  density  of  8  nodes  per 
meter  before  refinement. 
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Time  In  Seconds 


Figure  33.  SAFE-T  PML  Thickness  Analysis  (Vertical  Displacements)  of 
two  PMLs  one  3m  (i.e.,  3  wavelengths)  in  depth  and  the  other  2m  (i.e., 
2  wavelengths)  in  depth.  What  is  remarkable  in  this  analysis  is  how  well 
a  PML  as  thin  as  2m  compares  to  those  that  are  3m  and  beyond.  Note 
that  the  areas  of  greatest  change  are  not  at  the  peaks  and  valleys  of  the 
Rayleigh  wave.  This  is  a  huge  savings  computationally,  and  it  is  more 
efficient.  Table  III  illustrates  this  fact. 
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Figure  34.  SAFE-T  PML  Thickness  Analysis  (Total  Nodal  Displacement 
Magnitudes).  As  expected,  the  areas  of  greatest  change  are  not  at  the 
peaks  of  the  total  Rayleigh  wave  disturbance.  This  makes  the  PML  es¬ 
pecially  useful  in  analyzing  Rayleigh  waves  because  the  error  incurred  by 
thinning  the  PML  region  accumulates  away  from  the  area  of  interest,  that 
being  the  magnitude  of  the  Rayleigh  wave. 
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PML  Depth  Table 

Depth 

Elements 

Dofs 

Equations 

Memory 

Time 

5  meters 

12800 

26082 

78246 

147.6MB 

846.17  s 

4  meters 

10368 

21170 

63510 

126.0MB 

664.16  s 

3  meters 

8192 

16770 

50310 

98.9MB 

505.81  s 

2  meters 

6272 

12882 

38646 

79.4MB 

385.14  s 

1  meters 

4608 

9506 

28518 

134.9MB 

246.24  s 

Table  II.  This  PML  Depth  Table  displays  the  resource  cost  of  varying  the 
PML  thickness. 


PML  Efficiency  Table 

Depth 

Elements 

Dofs 

Equations 

Memory 

Time 

5  meters 

12800 

26082 

78246 

147.6MB 

846.17  s 

4  meters 

19.0% 

18.8% 

18.8% 

14.6% 

21.5% 

3  meters 

20.9% 

20.7% 

20.8% 

21.5% 

23.8% 

2  meters 

23.4% 

23.2% 

23.2% 

19.7% 

23.8% 

1  meters 

26.5% 

26.2% 

26.2% 

-69.9% 

36.1% 

Table  III.  This  table  calculates  the  benefit  of  reducing  PML  thickness. 
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IV.  MODEL  END-FIRE  ARRAY  OF 
SOURCE  THUMPERS 

A.  INTRODUCTION 

A  seismic  sonar  array  is  a  set  of  n  ground  source  elements  distributed  over 
an  area  of  the  Earth’s  surface  at  a  spacing  that  is  selected  to  allow  phasing  of  the 
excitation  of  individual  elements  to  constructively  or  destructively  contribute  to  a 
particular  source  radiation  pattern  [Ref.  46].  Modeling  of  a  linear  end-fire  sonar 
array  follows  the  principles  used  in  classical  antenna  configuration  theory.  The  total 
field  of  an  array  is  a  vector  superposition  of  the  fields  radiated  by  evenly  spaced  indi¬ 
vidual  elements.  Usually  array  elements  are  identical.  Directivity  can  be  achieved  by 
tuning  the  array  based  upon  its  geometric  configuration,  distance  between  elements, 
amplitude  and  phase  excitations,  and  the  radiating  patterns  of  the  individual  array  el¬ 
ements.  Research  on  seismic  SONAR  was  initiated  at  the  Naval  Postgraduate  School 
(NPS)  by  Dr.  Thomas  Muir  while  on  leave  from  the  Applied  Research  Laboratory  of 
the  LIniversity  of  Texas  at  Austin  (ARLTIT)  in  the  early  to  mid  1990s  with  the  goal 
of  determining  whether  buried  mines  could  be  detected  in  sand.  He  continued  this 
work  at  the  Naval  Postgraduate  School  when  he  became  chair  of  the  Aline  Warfare 
Department  in  the  late  1990s.  [Ref.  47] 

Prior  research  began  with  Lieutenant(USN)  William  Stewart.  He  was  first  to 
conduct  research  related  to  Seismo-Acoustic  SONAR  in  1995  [Ref.  48].  He  mounted  a 
plunger-type  source  using  a  loudspeaker  above  the  ground,  and  tested  the  transmitted 
signal  over  a  wide  range  of  frequencies.  His  tests  were  conducted  in  an  above  ground 
swimming  pool  filled  with  sand.  He  was  able  to  show  that  his  source  could  generate  a 
suitable  seismic  signal,  but  the  tank  was  too  small  for  any  echo  ranging  experiments. 

In  1998,  Lieutenant (USN)  Frederick  Gaghan  [Ref.  49]  focussed  his  research  on 
the  development  of  a  discrete-mode  excitation  source  that  consisted  of  two  inertial 
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mass  shakers  mounted  on  an  aluminum  framework.  Each  inertial  mass  shaker  was 
mounted  to  point  downward  at  a  45°  angle  in  an  effort  to  better  excite  elliptical 
Rayleigh  surface  waves.  While  promising  as  an  idea,  his  method  needed  a  more 
efficient  Rayleigh  wave  source.  Lieutenant (USN)  Sean  Fitzpatrick  [Ref.  50]  continued 
Gaghan’s  work  by  improving  on  the  source  using  a  linear  magnetic  force  actuator. 
With  a  two  element  seismometer  array,  he  was  able  to  locate  71-291  kg  targets  at 
ranges  up  to  5  meters  away.  Later  that  same  year,  student  Major (USMC)  Patrick 
Hall  [Ref.  51]  measured  the  reflectivity  of  targets  as  a  function  of  their  mass  load, 
and  found  that  target  reflected  signal  strength  was  proportional  to  target  mass. 

Captain(USA)  Kraig  Sheetz  [Ref.  52]  continued  seismic  SONAR  work  in  2000 
by  developing  a  receiver  that  was  capable  of  detecting  specific  objects  such  as  an  M-19, 
20  lb,  anti-tank  mine.  Lieutenant  (USN)  Scott  McClelland  [Ref.  53]  followed  Sheetz 
in  2002  by  mounting  two  inertial  shakers  onto  a  manually-pushed  rolling  cylinder. 
His  source  experiments  resulted  in  the  successful  detection  of  a  1000-lb  bomb  at  5 
meters.  Unfortunately,  the  roller  could  only  take  measurements  when  the  shakers 
were  directly  aligned  with  the  ground,  and  thus  it  proved  less  than  ideal. 

More  recently,  in  2003,  Lieutenant  (USN)  Douglas  MacLean  [Ref.  47]  intro¬ 
duced  a  small  tracked  vehicle  with  dual  inertial  mass  shakers  mounted  on  top  as  a 
mobile  source.  It  excited  Rayleigh  waves,  but  additionally  generated  unwanted  P- 
waves  that  destructively  interfered  with  signal  reception  of  surface  pulses,  thereby 
making  the  apparatus  incapable  of  finding  targets.  To  mitigate  the  destructive  in¬ 
terference  of  the  P-waves,  Lieutenant  (USN)  Steven  E.  Rumph  [Ref.  9]  developed  a 
four-element  end-fire  array  as  a  seismo-acoustic  SONAR  capable  of  being  spaced  and 
timed  in  such  a  way  as  to  constructively  interfere  Rayleigh  surface  waves  while  simul¬ 
taneously  destructively  interfering  unwanted  P-waves  and  body  waves.  Testing  on  a 
local  beach  yielded  3.5  meter  beam  patterns  with  approximately  15  db  suppression 
to  the  rear  of  the  array  relative  to  its  forward  direction. 
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Figure  35.  SAFE-T  element  by  element  deformation  for  medium  properties 
which  exhibits  a  compressional  wave  speed  of  1600  meters  per  second,  a 
shear  wave  speed  of  100  meters  per  second,  and  a  Rayleigh  wave  speed  of 
95  meters  per  second. 


Surface  Displacement 

of  a  node  9.5  meters  away 


Figure  36.  SAFE-T  vertical  displacement  results  with  no  PML  present. 
Time  is  represented  in  At  time-steps  which  are  0.0005  seconds  apart.  The 
medium  is  sand.  It  exhibits  a  compressional  wave  speed  of  1600  meters 
per  second,  a  shear  wave  speed  of  100  meters  per  second,  and  a  Rayleigh 
wave  speed  of  95  meters  per  second. 
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Surface  Displacement 

of  a  node  9.5  meters  away 


Figure  37.  SAFE-T  vertical  displacement  results  with  PML  present.  Notice 
the  smooth  transition  to  rest  on  the  surface  after  the  wave  passes.  Time  is 
represented  in  At  time-steps  which  are  0.0005  seconds  apart.  The  medium 
is  the  same  as  in  figure  36  which  exhibits  a  compressional  wave  speed  of 
1600  meters  per  second,  a  shear  wave  speed  of  100  meters  per  second,  and 
a  Rayleigh  wave  speed  of  95  meters  per  second. 

B.  SAFE-T  RESULTS:  OPTIMIZING  THE  AMPLITUDE 
OF  THE  SURFACE  WAVE 

Numerical  results  for  SAFE-T  are  presented  for  a  four-element  array  of  down¬ 
ward  (—z)  source  thumpers  on  a  half-space.  Figure  35  is  an  element  by  element 
deformation  graph  postprocessed  by  Altair  Engineering  Inc.’s  Hyperview  v7.0  using 
SAFE-T ’s  results  for  a  medium  with  properties  that  induces  a  compressional  wave 
speed  of  1600  meters  per  second,  a  shear  wave  speed  of  100  meters  per  second,  and  a 
Rayleigh  wave  speed  of  95  meters  per  second  inside  the  computational  domain.  Figure 
36  is  a  slice  of  the  numerical  domain.  Lateral  edges  have  rigid  Dirichlet  boundaries, 
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and  the  slice  is  taken  in  the  x  —  z  plane  of  the  domain.  In  the  first  instance,  there  is 
no  PML  present  to  absorb  outgoing  waves,  thus  unwanted  body  waves  bounce  back 
and  forth  between  the  rigid  boundaries.  The  material  properties  of  the  domain  are 
those  closely  related  to  sand,  namely,  Poisson’s  ration  v  =  0.498,  Young’s  Modulus 
E  =  8.06E07  Pa,  and  density  p  =  2690.00  kg/m3.  The  mesh  is  very  dense  (8  el¬ 
ements  per  meter)  in  order  to  provide  enough  nodes  to  minimize  dispersion  of  the 
source  pulse  using  At  time  steps.  Displacements  are  stably  computed  with  fourth 
order  accuracy  using  equation  11.110,  0.0005  seconds  apart.  The  medium  exhibits 
a  compressional  wave  speed  of  1600  meters  per  second,  a  shear  wave  speed  of  100 
meters  per  second,  and  a  Rayleigh  wave  speed  of  95  meters  per  second.  An  arbitrary 
Gaussian  source  is  used  to  initiate  a  four-element  end-fire  array  demonstrating  the 
effects  of  body  waves  in  a  medium  without  an  absorbing  layer.  Maximum  radiation 
occurs  in  the  direction  along  the  line  of  the  array  designated  as  the  positive  x  direction 
equivalent  to  0°  for  all  calculations.  SAFE-T  demonstrates  its  ability  to  effectively 
absorb  unwanted  body  waves  from  the  surface  of  the  computational  domain  in  figure 
37.  The  attenuation  component  of  the  damping  function,  equation  III.  16,  is  chosen  to 
be  linear  in  the  PML.  The  vertical  displacement  results  show  a  smooth  transition  to 
rest  on  the  surface  after  the  wave  passes.  Further,  because  of  the  speed  and  relatively 
small  amplitude  of  the  compressional  wave,  the  only  waves  clearly  visible  in  the  graph 
are  the  shear  and  Rayleigh  waves. 

Figure  38  shows  the  response  of  a  node  which  models  a  particles  motion  on 
the  free  surface  of  the  half-space  5  meters  from  the  end  of  the  end-fire  array.  It  is 
directly  in  the  path  of  maximum  radiation,  i.e.,  when  6  =  0°  or  along  the  axis  of 
the  array.  The  total  field  of  the  four-element  array  is  a  vector  superposition  of  the 
disturbances  generated  by  the  individual  element  thumpers.  In  order  to  provide  a 
more  directive  pattern,  it  is  necessary  to  have  the  partial  fields  (generated  by  the 
individual  thumpers)  interfere  constructively  in  the  direction  of  maximum  excitation 
and  interfere  destructively  in  the  remaining  wave  propagating  space. 


75 


Time  Delay  Comparisons 


5  meters  and  0  degrees  from  Array 


Figure  38.  This  graph  shows  the  response  of  a  node  located  on  the  free 
surface  and  5  meters  away  from  the  End-fire  array.  It  is  directly  in  the 
path  of  maximum  radiation,  i.e.,  when  6  =  0°  along  the  axis  of  the  array. 

a.  Time  Delay  and  Optimized  End-fire  Array 


Of  the  five  generally  excepted  methods  used  to  control  array  patterns, 
i.e.,  geometrical  configuration,  displacement  between  thumpers,  excitation  amplitude 
of  each  thumper,  excitation  phase  of  each  thumper,  relative  pattern  of  each  thumper, 
the  methods  found  most  effective  in  this  FE  model  were  space  between  thumpers 
and  phase.  The  phase  is  controlled  through  the  use  of  time  delay.  Time  delay  is 
accomplished  through  the  relation 


Time  Delay 


Distance  Between  Elements 
Wave  Speed 


(IV.  1) 


and  provides  an  effective  means  to  steer  through  interference  an  array  with  a  finite 
number  of  seismic  elements. 
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Time  Delay  Table 

Array 

Distance  Apart 

Wave  Speed 

Time  Delay 

Characteristics 

in  meters 

in  meters  per  second 

in  seconds 

4el  Linear 

0.25 

90  m/s 

0.002778  s 

4el  Linear 

0.25 

95  m/s 

0.002631  s 

4el  Linear 

0.25 

100  m/s 

0.002500  s 

4el  Linear 

0.25 

120  m/s 

0.002083  s 

4el  Linear 

0.25 

130  m/s 

0.001923  s 

4el  Linear 

0.25 

140  m/s 

0.001785  s 

Table  IV.  Time  delay  conversions  for  various  wave  speeds  allow  the  end-fire 
array  to  be  optimized  for  maximum  radiation  along  the  axis  of  propaga¬ 
tion. 


Table  IV  gives  time  delay  conversions  for  select  wave  speeds.  Figure 
39  shows  the  effects  of  destructive  interference  due  to  time  delay.  This  destructive 
interference  occurs  when  9  =  180°  along  the  axis  of  the  array  which  corresponds 
to  minimal  surface  wave  radiation.  There  is  a  marked  difference  in  the  amount  of 
destructive  interference  depending  on  the  time  delay  used.  In  this  model,  among  the 
three  delays  tested,  a  time  delay  of  0.002778  seconds  provides  the  most  destructive 
interference  behind  the  array. 

The  reaction  of  particles  to  body  waves  traveling  underneath  the  end- 
fire  array  is  of  primary  concern.  Wood’s(1968)  [Ref.  30]  experiments  show  that  there 
is  a  considerable  amount  of  energy  traveling  down  and  away  from  the  surface.  In 
order  to  optimize  the  energy  steered  by  the  array  0°  on  the  surface  and  along  the 
positive  axis  of  the  end-fire  array,  energy  traveling  down  must  be  minimized. 

Figure  40  depicts  the  nodal  displacement  magnitude  response  of  a  par¬ 
ticle  5  meters  underneath  the  end-fire  array  and  above  the  PML  in  a  non  layered 
media  .  The  largest  suppression  of  energy  occurs  when  a  0.002778  second  time  delay 
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Time  Delay  Comparisons 


5  meters  from  Array 


Figure  39.  This  graph  shows  the  response  of  a  node  located  on  the  free 
surface  and  5  meters  away  from  the  End-fire  array.  It  is  directly  in  the 
path  of  minimum  radiation,  i.e.,  when  9  =  180°  along  the  axis  of  the  array. 
The  effects  of  time  delay  clearly  shows  a  dramatic  reduction  in  the  amount 
of  energy  traveling  opposite  the  direction  of  steering. 


is  applied.  This  corresponds  to  a  wave  speed  of  90  meters  per  second.  The  other 
two  time  delays,  0.0025  seconds  (100  meters  per  second)  and  0.002083  seconds  (120 
meters  per  second),  also  show  a  suppression  of  body  waves  when  compared  to  the  non 
time-delayed  wave  strength.  Figure  41  shows  the  effects  of  constructive  interference  as 
the  timing  of  the  excitation  of  individual  elements  in  the  array  contribute  to  boosting 
the  wave’s  energy.  The  propagating  strength  of  the  wave  traveling  at  0°  and  along  the 
positive  axis  of  maximum  radiation  is  higher  than  the  non  time-delayed  wave.  Figure 
42  combines  into  one  graph  minimal  surface  nodal  magnitude  ( 6  =  180°),  maximal 
surface  nodal  magnitude  ( 6  =  0°),  and  the  nodal  displacement  effects  of  downward 
body  waves.  Figure  43  uses  the  derivative  of  the  Gaussian  wave  form  (figure  12)  as 
an  input  source  to  the  end-fire  array.  This  source  works  very  well  mathematically 


78 


Nodal  Displacement  Magnitude 


Time  Delay  Comparisons 


5  meters  from  Array 


Figure  40.  This  graph  shows  the  response  of  a  node  located  5  meters  under 
an  end-fire  array.  It  is  composed  of  mostly  shear  and  compressional  wave 
components.  Though  slight,  there  is  a  decrease  in  energy  propagating 
under  the  array  for  different  time  delays. 

because  it  generates  a  better  Rayleigh  wave.  By  specifically  tuning  the  source  to  a 
particular  Rayleigh  wave  speed  (see  Appendix  C),  the  derivative  of  the  Gaussian  best 
takes  advantage  of  the  PML  used  to  truncated  the  computational  domain.  Figure 
44  are  the  results  of  the  optimal  time  delay  for  achieving  the  highest  gain  at  0°  and 
along  the  positive  axis  of  the  array.  The  delay  is  0.002631  seconds  which  corresponds 
to  a  wave  speed  of  95  meters  per  second.  It  corresponds  to  the  minimum  radiation 
at  180°  and  under  the  array. 
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Time  Delay  Comparisons 


5  meters  from  Array 


Figure  41.  Time  delay  comparisons  taken  at  0  degrees  and  5  meters  away. 
It  includes  the  effects  of  time  delay  by  showing  the  propagating  energy 
when  no  time  delay  is  present. 


Time  Delay  Comparisons 

5  meters  from  Array 
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Figure  42.  This  graph  simultaneously  displays  time  delay  results  for  surface 
front,  surface  rear,  and  sub  surface  wave  propagation  nodal  displacement 
magnitudes. 
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Time  Delay  Comparisons 

3  meters  away 
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Figure  43.  The  use  of  the  derivative  of  the  Gaussian  pulse  as  a  time  delayed 
input  source  to  the  end-fire  array  takes  best  advantage  of  the  PML. 
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Figure  44.  The  optimal  time  delay  for  achieving  the  highest  gain  at  0° 
and  along  the  positive  axis  of  the  array  occurs  at  0.002631  seconds  which 
corresponds  to  a  wave  speed  of  95  meters  per  second.  The  times  for 
other  end-fire  time  delays  for  a  material  with  properties,  A  =  683.26P07  Pa, 
/i  =  2.69P07  Pa,  v  =  0.498,  Young’s  Modulus  E  =  8.06P07  Pa,  and  density 
p  =  2690.00  kg/m3  are  analyzed  and  compared. 
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V. 


FINDINGS  AND  CONCLUSIONS 


A  solid  time  dependent  perfectly  matched  layer  has  been  developed  to  absorb 
propagating  waves  which  result  from  a  seismic  event.  This  dissertation  presents  the 
major  steps  involved  in  building  a  transient  PML  model  for  an  isotropic,  homogeneous 
media  applied  to  truncating  the  computational  domain  where  elements  of  an  end-fire 
array  are  excited.  The  following  summarizes  the  findings  and  conclusions  of  this 
dissertation. 

•  Determined  method  to  find  suitable  analytic  benchmarks  for  seismic  wave 
analysis. 

In  order  to  measure  the  accuracy  of  any  finite  element  code,  suitable  bench¬ 
marks  have  to  be  analytically  computed.  For  the  seismic  pulse  on  a  half-space,  an¬ 
alytic  computations  involved  an  instantaneous  pressure  that  produces  a  singularity 
that  propagates  as  the  surface  wave.  Numerically,  this  presents  a  significant  chal¬ 
lenge.  How  does  one  determine  the  magnitude?  As  the  singular  point  is  approached, 
the  amplitude  of  the  disturbance  approaches  infinity.  As  a  result  of  this  dissertation, 
finite  analytic  solutions  exist  that  allow  numerical  methods  to  be  verified  for  accuracy 
and  efficiency. 

•  Development  of  a  new  strain-stress  equation  which  included  both  the  per¬ 
fectly  matched  media  and  the  computational  domain. 

This  dissertation  presented  the  first  known  three  dimensional,  vector- valued, 
time  dependent  stress-strain  relation  from  which  the  strain  and  stress  of  a  three  di¬ 
mensional  solid  system  within  a  perfectly  matched  medium  was  calculated.  This 
strain  equation  gives  the  damping  media  inhomogeneous  clastic  properties  that  at¬ 
tenuated  propagating  waves  in  the  PML  region.  Because  the  damping  properties  are 
dependent  upon  the  location  of  the  wave,  all  effects  of  the  attenuation  vanish  inside 
the  computational  domain. 

•  Development  of  3D  time  dependent  perfectly  matched  layer. 
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A  new  unsplit  3D  time  dependent  perfectly  matched  layer  was  derived  and  con¬ 
verted  into  its  weak  (Galerkin)  form  for  implementation  into  SAFE-T  (Solid  Adaptive 
Finite  Element  Transient),  developed  by  the  author  for  seismic  wave  analysis.  This 
dissertation  demonstrated  SAFE-T’s  ability  to  accurately  model  seismic  phenomena 
using  perfectly  matched  layers  to  absorb  unwanted  reflected  surface  and  body  waves. 

•  Determined  that  linearity  of  damping  function  did  not  contribute  greatly  to 
damping  region  properties. 

The  damping  functions  used  within  the  PML  region  determine  the  speed  at 
which  propagating  waves  are  attenuated.  A  comparison  of  the  effects  of  using  linear, 
quadratic,  and  cubic  damping  functions  showed  that  each  provided  excellent  absorp¬ 
tion.  The  analysis  did  not  reveal  a  significant  difference  in  the  rate  of  attenuation. 

•  Determined  that  the  damping  amplitude  when  made  too  great  has  the  effect 
of  stiffening  the  PML/computational  domain  interface  causing  reflections. 

The  damping  constant  stiffened  the  interface  by  essentially  making  what  should 
be  a  gradual  increase  in  attenuation  a  more  abrupt  change  thereby  causing  an  un¬ 
wanted  reflection.  The  analysis  of  this  dissertation  demonstrates  the  need  to  choose 
a  damping  constant  that  provides  maximum  amount  of  attenuation  with  the  least 
amount  of  stiffness  at  an  interface. 

•  Determined  that  reflections  cost  considerable  computer  resources  when  un¬ 
damped. 

Table  If  and  table  111  demonstrate  an  unintended  consequences  of  not  damping 
surface  and  body  waves,  namely, CPU  memory.  This  is  not  intuitive,  but  since  every 
motion  in  the  FE  model  need  be  calculated,  it  stands  to  reason  that  unwanted  reflec¬ 
tions  would  expend  computer  resources.  In  fact,  the  analysis  reveals  that  although  a  1 
meter  PML  requires  less  time  (a  mere  246.24  seconds  to  compute-compared  to  846.17 
seconds  for  the  5  meter  PML),  it  is  a  disaster  for  computer  memory.  It  consumes  an 
inordinate  amount  of  computer  memory.  It  requires  91.4%  of  the  memory  that  would 
be  needed  for  a  model  with  a  PML  5  meters  deep.  The  dissertation  demonstrates 
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further  the  need  to  attenuate  unwanted  reflections  within  a  FE  model. 

•  Determined  the  optimal  spacing  and  timing  for  maximal  Rayleigh  displace¬ 
ment  magnitude  and  minimal  body  wave  magnitudes  given  known  material  properties. 

SAFE-T  calculated  the  optimal  time  delay  and  space  between  elements  for 
achieving  the  highest  Rayleigh  surface  wave  gain  along  the  axis  of  the  end-fire  array. 
This  positive  axis  equates  to  0°  and  corresponds  to  the  end  of  the  array  that  produces 
the  largest  Rayleigh  wave.  Concurrently,  the  analysis  was  consistent  with  array  the¬ 
ory,  i.e.,  as  the  magnitude  of  the  Rayleigh  wave  increased,  SAFE-T  clearly  showed 
that  the  amplitude  of  the  body  waves  decreased. 

As  a  result  of  the  developments  involved  in  this  investigation,  several  future 
research  opportunities  exist.  Those  efforts  should  include,  but  not  limit  themselves 
to  the  following: 

•  Investigating  methods  to  make  the  PML  dynamic,  i.e.,  include  logic  into  the 
mc.ff  (MCOEFF)  FORTRAN  routine  to  sense  wave  motion  and  calculate  wave  speed 
for  the  damping  function. 

•  Place  obstacles  into  the  computational  domain  and  perform  scattering  cal¬ 
culations  and  source  level  estimations  on  a  variety  of  array  configurations. 

•  Analysis  of  non- homogeneous /anisotropic  materials. 

•  Analysis  of  non-linear  wave  phenomenon  such  as  shock  waves.  The  PML 
can  be  tuned  to  attenuate  non-linear  waves  as  well. 

•  Conduct  a  time  dependent  analysis  of  an  infinite  waveguide  using  the  tran¬ 
sient  PML  as  an  infinite  boundary. 

•  Examine  the  usefulness  of  method  in  non-destructive  testing  of  elastic  mem¬ 
bers  of  mechanical  devices  such  as  aircraft  wings,  nuclear  cooling  pipes,  and  ship  hull 
analysis. 
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APPENDIX.  A  (STRESS  TERMS) 

1.  UNKNOWN  STRESS  TERMS 


ja  (as®  f  +  f,s®  [f.<3+1  +  Fj<t'] )  Vl<man  (.1) 

a.  Main  Diagonal 

For  the  main  diagonal  the  terms  will  be  as  follows: 

TiV Va  =  S£  [(A  +  2/x)  SfjF,  rJJ1  +  AS|F2  r^1  +  AS3s3F3  ^J1]  (.2) 

^2,2  =  [AS^F,  <+!  +  (A  +  2/x)  Sf2F2  r^1  +  ASf3F3  r^1]  (.3) 

t?3+1V3,3  =  Si  [AS^F,  r";!1  +  AS|F2  r";+'  +  (A  +  2/x)  Sf3F3  R^1]  (.4) 

b.  Cross  Terms 

The  cross  terms  will  be  as  follows: 

rJ+%,2  =  S&,  (S®F,  «;«  +  SgF,  «JT)  (.5) 

t?s+\3  =  S&  (sf3F j  <J‘  +  S®  F3  ujj1)  (.6) 

t2”,+Vi  =  Sfo  (S®  F2  1  +  S®  F!  lift1)  (.7) 

t2"3+ Vs  =  S&,  (S|F2  +  S|F;)  «g‘)  (.8) 

TjYV.i  =  St>  (s®  F3  +  Sf3F3  »;«)  (.9) 

r32+V,2  =  S&,  (S|F3  <+‘  +  S|F2  «g')  (.10) 


2.  INITIAL  AND  KNOWN  STRESS  TERMS 


SfmXij  +  Ajm A t  J2  4  +  A2.„A«2  X)((n  +  1)  -  Oiy 


a.  Main  Diagonal 


For  the  main  diagonal  the  terms  will  be  as  follows: 

n  n 

Vi,i  =  -S^Xn  -  An  At  53  T11  -  An  A  t2  53((n  +  1)  -  0rn 

i=i  i=i 


(,ii) 


(.12) 
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(.13) 


1)2,2  —  —  $22X22  ~  A-22^-t  '^2  T22  —  A-22  At2  YX(n  +  1)  —  l) 


T. 


22 


1=1 


1=1 


v3,3  —  —  S33X33  —  A33  At  E  r33  -  a33a  t2  J2((n  + 1)  -  0 


T- 


33 


1=1 


1=1 


b.  Cross  Terms 

The  cross  terms  will  be  as  follows: 


vi,2  —  —S22X12  —  A22A t  Y  ri2  ~  A-22  At2  Y((n  +  1)  —  0T12 


1=1 

n 


1=1 

n 


Vl,3  =  -S33X13  -  A33  At  Z  Ti3  -  A33A t2  +  1)  -  l) 


T- 


13 


1=1 


1=1 


V2,i  =  -SnX2i  -  An  At  Y  r2i  -  An  At2  ^((n  +  1)  -  l) 


T , 


21 


1=1 


1=1 


v2,3  ~  “833X23  —  A33  At  Yz  r23  -  A33A t2  J2((n  +  !)  -  0 


T , 


23 


1=1 


1=1 


v3,i  =  -SnX3i  -  An  At  ^  rl31  -  An  At2  Y((n  +  1)  -  l) 


T. 


31 


1=1 


1=1 


v3,2  ~  822X32  ~  A22  At  Y  r32  -  a22a  t2  J2((n  + !)  -  0 


T- 


32 


1=1 


1=1 


(.14) 


(.15) 

(.16) 

(.17) 

(.18) 

(.19) 

(.20) 


INITIAL  AND  KNOWN  BOUNDARY  TERMS 


'r, 


i?(x,y,z)vids 


a.  Three  Components 


(.21) 


(.22) 

(.23) 

(.24) 
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v  1  =  ti(x, y,  z) 
V2  =  i%(x,y,  z) 
n  =  il(x,y,z) 


APPENDIX.  B  (MCOEFF  TERMS) 


1.  MCOEFF  TERMS 

Subroutine  MCOEFF  is  called  for  every  integration  point  and  evaluates  all  of 
the  components  of  the  coefficient  matrices  for  a  single  element  or  a  batch  of  elements. 
The  AtJ  PHLEX  coefficients  are  identified  from  the  final  Galerkin  (weak)  form  of  the 
PDE. 


Set  1 


Set  2 


a.  LHS  Coefficients 

Body  forces  induced  as  time  marches. 


(M1  Vi)A00(l,l)  —  Ki  +  fk  +  ^ 

(.25) 

(u2  v2)A00(2,  2)  —  K\  T  fk  +  ^ 

(.26) 

(u3  p3Moo(3,  3)  —  +  fk  +  ^ 

(.27) 

Vi)-4„(1, 1)  =  (A  +  2/x)  SftSftF! 

(.28) 

nl+2'vl,2)A22(l,l)  =  flS-*2S?2F1 

(.29) 

5+1t.1,i)^12(l,2)  =  ASf1Sf1F1 

(.30) 

5!1i-i.2)Ai(l,2)=^2Sf1F2 

(.31) 

51«u)A33(l,l)=/2S4sf3F1 

(.32) 

JJV.i)Ali3(l,3)  =  ASf1Sf1F1 

(.33) 

5,l1i-i.3)Ai(l,3)=/.SiSf1F3 

(.34) 

^v2,1)A11(2,2)  =  fiSt1SgF2 

(.35) 

$%,2)  A22(2, 2)  =  (A  +  2/i )  S£S?2F2 

(.36) 

-+1u2>i)2412(2,l)  =  ^S(\Sf2F1 

(.37) 
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«+%,2)A2  i(2,  1)  =  \S*SgF2  (.38) 

K+1V2)3)^33(2,2)  =  //S4S^F2  (.39) 

«+1U2)2)^23(2, 3)  =  AS£S?2F2  (.40) 

Kt1^2,3)-432(2, 3)  =  /xS4s|F3  (.41) 

Set  3 

KJ'os.iMiaP, 1)  =  JiSf1Sf3F1  (.42) 

K+VMs.P,  1)  =  AS^S|F3  (.43) 

(“S'l’s, 2)4123(3, 2)  =  (iS4s|F2  (.44) 

KJ1f3.3)4l32(3, 2)  =  AS^Sf3F3  (.45) 

KI'fs.iMnP,  3)  =  AiSfjSfjF,  (.46) 

K««3,2)4l22(3, 3)  =  /(S.22Sf2F3  (.47) 


(“S'l’s, 3)4133(3, 3)  =  (A  +  2/3)  S4S|F3  (.48) 

b.  RHS  Coefficients 

(vi)f(l)  =  Kl<  +  «2^  +  -  flAti^U1!  (.49) 

r)vn  f)2vn  n 

(u2)/(2)  =  Ki<  +  K2^  +  -  flAtJ2U2  (-50) 

r)iin  fPlln  ^ 

(u3)/(3)  =Ki«3  +  «2^-  +  k3-^  -  flAt^Us  (-51) 

n  n 

(vi,i)fx(l)  =  ~S?lXu  -  An  A  t  X  rh  -  An  At2  X((n  +  !)  -  0rn  (-52) 

1=1  1=1 

n  n 

(ui,2)/®(2)  =  -S4xi2  -  A22At  X  ri2  -  A22At2  X((™  +  1)  -  0ri2  (-53) 

Z=1  Z=1 

n  n 

(ui,3)/®(3)  =  -S33Xi3  -  A33At  X  ris  -  A33At2  X((w  +  x)  -  0rf3  (-54) 

1=1  1=1 

n  n 

(v2,l)fy(l)  =  -SfiX21  -  An  At  X  r2i  -  AnAt2  X((w  +  !)  -  0r21  (-55) 
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(V2,2)fy(^)  ~ 

(V2,3)fy{3)  = 

Kl)/2(1)  = 

(^3,2) /2  (2)  = 

(V3,3)fz(  3)  = 


-S22X22  -  A22At  X  r22  -  A22At2  X((™  +  1)  -  /)r22  (.56) 


«=i 


z=i 


S33X23  -  A33At  X  r23  -  A33At2  X((™  +  1)  -  /)r^3  (.57) 

/=i  z=i 


-SnXsi  -  An  At  X  r3i  “  An  At2  X((n  +  1)  -  l) 
1=1  1=1 

n  n 

-$4x32  —  a22a  t  X  r32 —  a22a t2  X((n  + 1)  ~  0 

1=1  1=1 

n  n 

■8.4x33  —  A33A t  X  r33  —  A33At^  X((n  +  1)  —  0 


31 


32 


33 


(.58) 

(.59) 

(.60) 


1=1 


1=1 
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APPENDIX.  C  (SOURCE  COMPUTATIONS) 


The  following  computer  code  found  in  this  appendix  presents  the  procedures 
used  to  construct  the  source  used  for  the  seismic  sonar  array,  ft  was  crafted  in 
such  a  way  as  to  produce  a  Rayleigh  wave  of  unit  wavelength  for  material  properties 
which  match  closely  with  that  of  sand.  It  is  written  and  evaluated  using  Wolfram’s 
Mathematica  5.2. 


F(W) 
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Source  Transforms2.nb 


1 


In [1]  :  = 
In  [21:  = 


Out [6]= 
Out  [7]  = 


Out [8]= 


Out  [9]  = 

In  [10]  r 
Out  [10]: 

In  [11]  :  - 
Out  [11  ]: 


Of f [FindMaximum: : "lstol",  General: : "spelll"] 

Amplitude  =  1 ; 

Peak  =  0.005500; 

Wid  =  0 . 0023583; 

f[t„]  :=  Amplitude  ®^<t-psak)''2/widA2 

f  [t] 
f ’  [t] 

ftplot  =  Plot  [f  [t]  ,  {t,  0,  2  Peak},  PlotRange ->  {  0 ,  Amplitude}] 
derftplot  =  Plot  [f  1  [t]  ,  {z,  0,  2  Peak}  ,  PlotRange  ->  {-500 . 0,  500.0}] 


e-179805.  (-0 . 0055+t) 2 

-359610.  e-179805-  (-0. 0055+t)2  (-0.0055  +  t) 


-  Graphics  - 


-  Graphics  - 

FindRootff  1  [t]  ==0,  {t,  0.001,  0.006}] 

{t  ->  0.0055} 

FindRoot  [f  ’  ’  [t]  ==  0,  {t,  0,  0.004}] 

{t  >4  0.00383243} 
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Source  Transforms2.nb 


2 


In [12]:= 
Out [12]= 

In  [13]  :  = 
Out [13]= 

In [14] := 
Out  [14]  = 

In [15] := 
Out  [15 ]  = 

In  [16]  :  = 


Out  [16]  = 
In  [17]:  = 
Out  [17]  = 
In  [18]  :  = 
In  [19]  :  = 
Out  [19]  = 


In  [20]  :  = 
Out [20]= 


anst  =  t  /  .  % [ [1] ] 

0.00383243 


f ' [anst] 

363 . 721 


Solve  [h  f  '  [anst]  ==  1,  h] 

{{h->  0.00274936}} 

ansh  =  h  /  .  %[  1] ] 

0.00274936 


derftplot  =  Plot  [ansh  f  1  [r]  ,  [z,  0,  2Peak},  PlotRange  -»  {  - 1 ,  1}] 


Graphics 


ansh  f  '  [r] 

-988.697  e~119805'  (-0-0055+r)2  (_o.0055  +  t) 

utrans[£d_]  :  =  FourierTransform[f  1  [t]  ,  t,  w] 
utrans [w] 

-143464. 

(o  .  e-8*67362*10-19  iw+0-  i0j2  +  @1  •  39039xl0-6  ( 1977 . 86+i  w)  2  (2.94319  X  10“23  +  (0  .  - 
ei.39039xicr6  (1977.86+iw)2  (  ( _  4 . 9  9  2  6  9  x  1 0  ~8  Erf  [-2.33219  -  0.00117915 
4. 99269  x  10~8  Erf  [2.33219  +  0.00117915ito]  Sign  [1977. 86  h 
(  (-4.99269X  10~8  +  0.  i)  -  (0.  +  2.52429x  10'11  i)  to) 

(Erf [-2.33219  -  (0.  +  0 . 0011 7915  i )  to]  + Erf [ 2 . 332 19  +  (0.  +  0. 
Sign  [2. 33219  +  (0.  +  0.00117915i)  6J  ]  2 ) 

Abs [utrans [1] ] 2  /  N 

2 . 78078  x  10~6 


+  5 . 04859  x  10-11  i)  to)  + 

i  to]  - 

-iu]2t 

00117915  i )  to]  ) 
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Source  Transforms2.nb 


3 


In  [21  ]  :  =  Plot [Abs [utrans [w] ]  ,  {(*) ,  0,  2000},  PlotRange  -»  Automatic,  PlotPoints  -»  100] 


0ut[21]=  -  Graphics  - 

In  [22  ]  :  =  FindMaximum  [Abs  [utrans  [a]  ,  {&>,  500}] 

Out  [22  ]  =  {0.606531,  {0)^  599.675}} 

In  [23]  :=  6J  =  599 . 6749964809572' 

0ut[23]=  599.675 

Looking  for  a  wavelength  of  1  meter.  Use  =  /  and  then  =  wavelength  *  cR 


In  [24]:=  - 

2  7T 

Out  [24]  =  95.4412 

Material  properties  for  sand... 

kg 


In  [25]:=  p  =  2690 


H  =  2.69*10 


v 


Out [25 ] = 


X =  683.26*10 


2690  kg 


s^  m 
kg 


Determine  Poisson's  Ratio  given  the  material  properties... 

A 

In  [28]:=  v  = 


2  (X  +  u) 

Out  [28  ]  =  0.498039 
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Source  Transforms2.nb 


4 


In  [29]  :  = 


:  =  cs  =  PowerExpandJ  — J  ] 


Out  [29]  = 


Computing  the  Rayleigh  wave  speed  using  Achenbach... 


.862  +  I.14v 

In[30]  :  =  cR  =  -  cs 

1  +  v 


Out [30]= 


95.4424m 
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APPENDIX.  D  (ARRAY  SUPERPOSITION 

CALCULATIONS) 


This  appendix  gives  the  Mathematica  code  used  to  examine  the  analytic  prop¬ 
erties  of  an  end-fire  array. 


Surface  Response  1  m  away 


Surface  Source 


Pekeris  1 955 

0.3  I - 1 - 1 - 1 - 1 - T - 1 - 1 - 1 - 1 - 

o.2  - . .  . . . ;• . .  . .  . - 

0.1- . .  . . .  . .  . .  . - 

c  i 

CD  ......... 

E  0-  . . . . .  . .  . - 

CD  | 

CO  \ 

q.  -o.i  - . . ; . . . 'A- 1 . . ; . . ; . - 

Q  \  : 

-0.2  - . .  . . . . .  . .  . - 

-o.3  - . ; . ■; . ; . ; . ;■■■■ . ; . ;■ . ; . ■; . - 

-0.4 - ' - ' - ' - ' - 1 - ' - ' - ' - 1 - 

0  0.2  0.4  0.6  0.8  1  1.2  1.4  1.6  1.8  2 

time  (seconds) 
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EndfireArray2.nb 


1 


Endfire  1 /4-lambda  Array  w  /  Gaussian 
Source 


In [1]  :  = 


Amplitude  =  1 ; 
Peakl  =  .0055; 
Widl =  .0023583; 
Peak2  =  .0080; 
Wid2 =  .0023583; 


Peak3  =  .0105; 
Wid3 =  .0023583; 


Peak4  =  .0130; 


Wid4 =  .0023583; 
f 1 [t_]  :  =  Amplitude  © 
f2  [t_]  :  =  Amplitude  © 
f 3 [t_]  :  =  Amplitude  © 
f4  [t_]  :=  Amplitude© 
ftplotl  =  Plot[fl[t] , 
ftplot2  =  Plot [f 2 [t] , 
ftplot3  =  Plot [f3[t] , 
ftplot4  =  Plot [f4 [t] , 


- (t-Peakl) A2/WidlA2 
- (t-Peak2) A2/Wid2A2 
- (t-Peak3) A2/Wid3A2 
-  (t-Peak4) A2/Wid4A2 

{t,  0,  5  Peakl},  P  lot  Range  -»  {0,  Amplitude}] 

{t,  0,  5  Peakl},  P  lot  Range  -»  {0,  Amplitude}] 

{t,  0,  5  Peakl},  P  lot  Range  -»  {0,  Amplitude}] 

{t,  0,  5  Peakl},  P  lot  Range  -»  {0,  Amplitude}] 


Out  [14]  =  -Graphics 
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Endfire Array!,  nb 


Out[15]=  -Graphics 


Out[16]=  -Graphics 


Out[17]=  -Graphics 


EndfireArrayl.nb 


3 


fl[t]  +f2[t]  +  f 3  [t]  +f4[t] 

In  [18]  :=  gl[t_]  :=  - 

4 

Plot[gl[t],  {t,  0,  0.025},  PlotRange ->  { 0 ,  Amplitude}] 


Out[19]=  -  Graphics  - 

In  [20]  :=  utrans  [&) _ ]  :  =  FourierTransf  orm[gl  [t]  ,  t,  0)] 

utrans [u] 

- i=  (0 . 0000181552  (g1- 39039x10-6  +  9 . 07761  x  10-6  e1'39039"10'6  (iW-W+i-)* 

4  V2tt  ' 

(1.  Erf  [-2.33219  -0.0011 7915iw]  +1.  Erf  [2.33219  +  0.0011 7915iw] 

Sign  [1977. 86  +  iw]2  +  2.10131xl0-8  e1-39039*'10-6  (2876.88+iu)2 

(2  .  +  (1 .  Erf  [-3 .39227  -  0 . 0011 7915  iw]  +1.  Erf  [3 .39227  +  0.00117915iw]) 
Sign  [2876 . 88  +  i  w]  2)  +  5 . 1393  x  10-12  e1'39039"10"6  <3775-91+i“>2 
(2  .  +  (1 .  Erf  [-4. 45236  -  0 . 0011 7915  iw]  +1.  Erf  [4. 45236  +  0.00117915iw]) 
Sign  [3775 . 91  +  i  u]2)  +  1 . 32805  x  10-16  e1'39039x10-6  <4^74 . 93+i a.) 2 
(2  .  +  (1 .  Erf  [-5.51245  -  0. 0011 7915  iw]  +1.  Erf  [5.51245  +  0.00117915iw]) 
Sign [46  74.93  +  i  w]2)  ) 

In  [21  ]  :  =  Abs  [utrans  [  .  01]  ]  2  /  N 

Out  [21  ]  =  1 . 08932  x  10-6 

In  [22  ]  :  =  Plot  [Abs  [utrans  [to]  ]  2 ,  {01 ,  0,  2000},  PlotRange  -»  Automatic,  PlotPoints  -»  100] 
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EndfireArrayl.nb 


4 


Endfire  1/4-lambda  Array  w/  Derivative  of 
Gaussian  Source 

In  [1  ]  :  =  Amplitude  =  1; 

Peakl  =  .0055; 

Widl =  .0023583; 

Peak2  =  .0080; 

Wid2 =  .0023583; 

Peak3  =  .0105; 

Wid3 =  .0023583; 

Peak4  =  .0130; 

Wid4 =  .0023583; 

fl[t_]  :=  Amplitude  ®-<t-peakl)''2/wid1''2 

f2  [ t _ ]  :=  Amplitude  ®-<t“Paak2)''2/wid2''2 

f  3  [ t _ ]  :=  Amplitude  ®-(t-paak3),'2/wid3''2 

f4[t_J  :=  Amplitude  ®-<t“Peak4)''2/wid4''2 

ftplotl  =  Plot  [f  1  '  [t]  ,  {t,  0,  5  Peakl},  PlotRange ->{- 500  Amplitude,  500  Amplitude}  ] 

ftplot2  =  Plot  [f2 '  [t]  ,  {t,  0,  5  Peakl},  PlotRange  ->{- 500  Amplitude,  500  Amplitude}  ] 

ftplot3  =  Plot  [f  3 '  [t]  ,  {t,  0,  5Peakl},  PlotRange ->{- 500  Amplitude,  500  Amplitude}  ] 

ftplot4  =  Plot  [f4 '  [t]  ,  {t,  0,  5  Peakl},  PlotRange  ->{- 500  Amplitude,  500  Amplitude}  ] 


0ut[14]=  -Graphics 
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EndfireArrayl.nb 


Out[15]=  -Graphics 


Out[16]=  -Graphics 


Out[17]=  -Graphics 
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APPENDIX.  E  (PROPHLEX  FORTRAN 

MODULE) 


c  —  mc.ff  -  Fri  Jul  26  14:29:13  CDT  2006 

c  Copyright:  Computational  Mechanics,  Co.,  Inc.  1992-1996 
c  MAJ  Anthony  N.  Johnson,  Naval  Postgraduate  School 
c  Monterey,  California  93940 


#include 

"PH-stdF.h" 

#include  ' 

'PH-parameters . h" 

#include  ' 

'PH-Fmacros .h" 

#include  ' 

'application-parameters .h" 

#include  ' 

c 

'application-macros .h" 

c 

subroutine  mc(el2get, 


& 

a00 , aOl , a02 , a03 , 

k 

al0,all ,al2,al3. 

k 

a20 , a21 , a22 , a23 , 

k 

a30 , a31 , a32 , a33 , 

k 

f ,fx,fy ,fz, 

k 

xyz,  xyznod,  ngnode, 

k 

dxdxi , dxidx , x j  ac , 

k 

numel ,ncomp ,nca,nsol 

k 

u , uxyz , 

c 

k 

xi ,eta,zeta,bcinf o) 

C*  * 

c*  Routine:  me  * 
c*  Purpose:  define  the  interior  integral  coefficients  of  the  * 
c*  variational  problem.  * 
c*  Variables:  * 

c*  -  * 

c*  I  el2get :  element  number  to  load  if  >  0  * 
c*  otherwise  batch  flag  if  <=  0  * 
c*  0  a00,...,a33:  coefficients  of  the  contribution  to  the  stiffness* 
c*  matrix  from  the  interior  integral  coded  as:  * 
c*  apq:  (p,q)  =  {0, 1,2,3}  where:  * 
c*  0:  signifies  the  shape  function  itself  * 
c*  1:  signifies  x-derivative  of  shape  function  * 
c*  2:  signifies  y-derivative  of  shape  function  * 
c*  3:  signifies  z-derivative  of  shape  function  * 
c*  p:  first  index  of  apq  signifies  test  f unc .  * 
c*  q:  second  index  of  apq  signifies  trial  func-* 
c*  tion  (or  the  solution  and  its  derivs . )  * 
c*  0  f,...,fz:  coefficients  of  the  contribution  to  the  load  * 
c*  vector  of  the  interior  integral  coded  as:  * 
c*  fp:  p  =  {  ,n,t,s}  where:  * 
c*  :  signifies  shape  function  itself  * 
c*  x:  signifies  x-derivative  of  shape  function  * 
c*  y:  signifies  y-derivative  of  shape  function  * 
c*  z:  signifies  z-derivative  of  shape  function  * 
c*  p:  signifies  a  function  multiplied  by  the  * 
c*  test  function  and  its  derivatives.  * 
c*  I  xyz:  integratin  point  coords  * 
c*  I  xyznod:  nodal  point  coordinates  for  each  element  * 
c*  I  ngnode:  number  of  nodes  per  element  * 
c*  I  dxdxi:  dx/dxi  * 
c*  I  dxidx:  dxi/dx  * 
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c* 

I 

xjac 

jacobians  for  each  element 

c* 

I 

numel : 

number  of  elements  in  the  batch 

c* 

I 

ncomp : 

number  of  solution  components 

c* 

I 

nca: 

number  of  active  solution  components 

c* 

I 

nsol : 

number  of  solutions 

c* 

I 

u , uxyz : 

solution  and  its  derivatives  at  the  intgr.  point 

c* 

I 

xi ,eta,zeta 

integration  point  coords  in  master  coords 

c* 

p 

bcinf o 

boundary  condition  info 

c* 

Q^e3|e^c3f;3|e^c3f;^e^c3|e3|e^c3|e^e3f;3|e^e3f;3|e^c^e^e^c3)c3|e^c3|c^e^c3|e^e3f;3|e^e3f;3|e^c3f;3|c^c3f;^e^c^e^e^c3|e3|e^c3|e^e3f;3|e^c^e3|e^c3f:3|e^c^e^e^c3|e3|e^c3|c^e^e3|e^e 

c 

#include  "PH-std.blk" 

#include  "PH-ftnwrk.blk" 

#include  "PH-wrkspcPtrs .blk" 
c 

REALTYPE 

Sc  a00,a01  , a02,a03. 

Sc  al0,all,al2,al3. 

Sc  a20 ,  a21 ,  a22 ,  a23 , 

Sc  a30 ,  a31 ,  a32 ,  a33 , 

&  f,fx,fy,fz, 

Sc  xyz,  xyznod, 

&  dxdxi,dxidx,xjac, 

Sc  u,uxyz,xi,eta,zeta 

c 

REALTYPE  xlambdat ,  xmut ,  xrhot ,  pi 
REALTYPE  beta,  gamma,  deltat 
REALTYPE  fprp,  fevn, 

Sc  pmlfx,  pmlfy,  pmlfz,  cs 

REALTYPE  fm,fc,fk,fl,kl,k2,k3, 

Sc  tel  ,tc2,tc3,tc4,tc5,tc6 

REALTYPE  sum_e,  sum_ui j ,  sum_u, 

Sc  sum_T,  sum2_T,  sum2_e, 

Sc  etn,  etnpl,  Ecur,  Tcur 

REALTYPE  Lpmlx,  Lpmly,  Lpmlz 
REALTYPE  mx,  my,  mz,  bx,  by,  bz 
REALTYPE  pmlstartx,  pmlstarty,  pmlstartz 
c 

INTTYPE  el2get,  ngnode,  numel,  nod 
INTTYPE  ncomp,  nca,  nsol,  lastepl 
INTTYPE  soltn,  soltnpl,  veltn,  acctn 
INTTYPE  k,  1,  lastep2 
c 

CPTRTYPE  elptr ,bcinfo 
c 

INTTYPE  ielstrt,  ielend,  curstep 
INTTYPE  icomp,  jcomp,  iel 
INTTYPE  debug,  pml_on 
INTTYPE  passnum,  iresid,  idiag 
c 
c 

dimension 

Sc  aOO  (numel ,  ncomp ,  ncomp)  , 

Sc  aOl  (numel ,  ncomp ,  ncomp)  , 

Sc  a02  (numel ,  ncomp ,  ncomp)  , 

Sc  a03 (numel, ncomp, ncomp)  , 

Sc  alO (numel, ncomp, ncomp)  , 

Sc  all  (numel , ncomp , ncomp)  , 

Sc  al2  (numel , ncomp , ncomp)  , 

Sc  al 3 (numel, ncomp, ncomp)  , 

Sc  a20  (numel ,  ncomp ,  ncomp)  , 

Sc  a21  (numel , ncomp, ncomp )  , 

Sc  a22 (numel, ncomp, ncomp)  , 

Sc  a23 (numel, ncomp, ncomp)  , 
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&  a30  (niimel ,  ncomp ,  ncomp)  , 

&  a31 (numel , ncomp , ncomp) , 

&  a32 (numel , ncomp , ncomp) , 

&  a33 (numel, ncomp, ncomp) 

dimension 

&  etn (MAXELEMBATCH, ncomp, ncomp) , 

&  etnpl (MAXELEMBATCH , ncomp , ncomp) , 

&  Ecur (MAXELEMBATCH, ncomp, ncomp) , 

&  Tcur (MAXELEMBATCH, ncomp, ncomp) , 

&  sum_e (MAXELEMBATCH , ncomp , ncomp) , 

&  sum2_e (MAXELEMBATCH, ncomp, ncomp) , 

&  sum_T (MAXELEMBATCH , ncomp , ncomp) , 

&  sum_ui j (MAXELEMBATCH , ncomp , ncomp) , 

&  sum2_T (MAXELEMBATCH , ncomp , ncomp) 

dimension 

&  u (numel , ncomp ,nsol) , 

&  sum_u (MAXELEMBATCH, ncomp, nsol ) , 

&  uxyz (numel, ncomp, nsol, 3) , 

&  f (numel , ncomp) , 

&  f prp (MAXELEMBATCH , ncomp) , 

&  f evn (MAXELEMBATCH, ncomp) , 

&  f x (MAXELEMBATCH , ncomp) , 

&  f y (MAXELEMBATCH , ncomp) , 

&  fz (MAXELEMBATCH, ncomp) 

dimension 

&  xyz (numel , 3) ,  xyznod (numel , ngnode , 3) , 

&  dxdxi(numel,3,3) , dxidx (numel ,3,3) ,xjac (numel) , 

&  bcinf o (numel ,7) 

dimension 

&  xlambdat (MAXELEMBATCH) , 

&  xmut (MAXELEMBATCH) , 

&  xrhot (MAXELEMBATCH) , 

&  fm (MAXELEMBATCH) , 

&  fc (MAXELEMBATCH) , 

&  fk (MAXELEMBATCH) , 

&  fl (MAXELEMBATCH) , 

&  kl (MAXELEMBATCH) , 

&  k2 (MAXELEMBATCH) , 

&  k3 ( MAXELEMBATCH ) , 

&  tel (MAXELEMBATCH) , 

&  tc2 (MAXELEMBATCH) , 

&  tc3 (MAXELEMBATCH) , 

&  tc4 (MAXELEMBATCH) , 

&  tc5 (MAXELEMBATCH) , 

&  tc6 (MAXELEMBATCH) , 

&  pmlfx (MAXELEMBATCH) , 

&  pmlfy (MAXELEMBATCH) , 

&  pmlfz (MAXELEMBATCH) , 

&  cs (MAXELEMBATCH) 

c 

pml_on  =  1 
c  debug  =  1 

c - 

c  load  the  time  step  parameter  data 
c 

call  GTAPPROP(beta,gamma) 

call  GTPARAM (  NONL_PARAMS ,  DTNONL,  deltat) 

call  GTPARAM (NONL_PARAMS,  ITSTEP,  curstep) 

c -  set  up  the  parallel  batch  stuff 

c 

call  PH_GET_BATCH_(  el2get,  numel,  ielstrt,  ielend  ) 
c 
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c -  get  residual  flag,  the  diagonal  preconditioner  flag,  and 

c  the  pass  number  flag 

c 


PH_GET_SWITCH_ (S_resid ,  iresid) 

PH_GET_SWITCH_ (S_diag ,  idiag) 

PH_GET_SWITCH_(S_passnum,  passnum) 
c 
c 

call  PH_UG_GET_SLOT_ (S0LUTI0N_TN ,  soltn) 
call  PH_UG_GET_SLOT_ (VELOCITY_TN ,  veltn) 
call  PH_UG_GET_SLOT_ (ACCELERATION_TN ,  acctn) 
c 

c  load  up  the  material  data  for  the  element  batch 
c 

do  3101  iel  =  ielstrt , ielend 
c 

elptr  =  PH_GET_ELEMENT_POINTER_FROM_BATCHID_(iel) 
call  gethooke  (  elptr,  xlambdat (iel) ,  xmut(iel),  xrhot(iel), 
&  pmlfx(iel),  pmlfy(iel) ,  pmlfz(iel)  ) 

c 

if  (xrhot(iel)  .eq.  0)  then 
cs (iel)=0 . OdO 

else 

cs (iel) =sqrt (xmut (iel) /xrhot (iel) ) 
endif 

c 

do  3111  nod  =  1,  ngnode 
c 

if  (abs (xyznod(iel ,nod, 1) )  .gt.  Lpmlx)  then 
Lpmlx  =  abs (xyznod(iel ,nod, 1) ) 
c  print  *,  ’Lpmlx’,  Lpmlx 

c  pmlstartx  =  Lpmlx 

endif 
c 

if  (abs(xyznod(iel,nod,2))  .gt.  Lpmly)  then 
Lpmly  =  abs (xyznod(iel ,nod,2) ) 
c  print  *,  ’Lpmly’,  Lpmly 

c  pmlstarty  =  Lpmly 

endif 
c 

if  (abs (xyznod(iel ,nod,3) )  .gt.  Lpmlz)  then 
Lpmlz  =  abs (xyznod(iel ,nod,3) ) 


c  print  *,  ’Lpmlz’,  Lpmlz 

c  pmlstartz  =  Lpmlz 

endif 
c 

3111  continue 

3101  continue 

c 

do  4001  iel  =  ielstrt , ielend 
c -  calculate  propagating  wave  functions 


c 

c 

fprp(iel,l)  =  abs(pmlfx(iel)*(xyznod(iel,2, 1) /Lpmlx) **2) 
fprp(iel,2)  =  abs (pmlf y (iel) * (xyznod(iel , 2 , 2) /Lpmly) **2) 
fprp(iel,3)  =  abs(pmlfz(iel)*(xyznod(iel,2,3)/Lpmlz)**2) 
c 

c -  calculate  evenescant  wave  functions 

c 

fevn(iel,l)  =  abs (pmlf x(iel)*(xyznod(iel, 2, 1) /Lpmlx) **2) 
fevn(iel,2)  =  abs(pmlfy(iel)*(xyznod(iel,2,2)/Lpmly)**2) 
fevn(iel,3)  =  abs(pmlfz(iel)*(xyznod(iel,2,3)/Lpmlz)**2) 
c 


no 


c 

4001  continue 
c 

c -  initialize  coefficients 

c  note:  do  not  initialize  coefficients  which  are  not  active 

c  ie.  if  al3  is  not  active  omit  the  initialization,  this 

c  will  overwrite  memory! ! ! ! 

c 


do  5001  iel  =  ielstrt , ielend 
do  5011  icomp  =  l,nca 
do  5021  jcomp  =  l,nca 

aOO (iel , icomp, j comp)  =  O.OdO 
all (iel , icomp, jcomp)  =  O.OdO 
al2 (iel , icomp, jcomp)  =  O.OdO 
al3 (iel , icomp, jcomp)  =  O.OdO 
a21 (iel , icomp, jcomp)  =  O.OdO 
a22 (iel , icomp, jcomp)  =  O.OdO 
a23 (iel , icomp, jcomp)  =  O.OdO 
a31 (iel , icomp, jcomp)  =  O.OdO 
a32 (iel , icomp, jcomp)  =  O.OdO 
a33 (iel , icomp, jcomp)  =  O.OdO 
Ecur( iel, icomp, jcomp)  =  O.OdO 
Tcur (iel , icomp, jcomp)  =  O.OdO 
5021  continue 
5011  continue 
5001  continue 
c 
c 

do  6101  iel  =  ielstrt , ielend 

c -  compute  necessary  coefficient  values 

c 

fm(iel)  =  xrhot(iel)  *  (l.OdO  +  fevn(iel,l)) 

&  *  (l.OdO  +  fevn(iel,2))  *  (l.OdO  +  fevn(iel,3)) 

fc(iel)  =  xrhot(iel)  *  cs(iel)  * 

&  (  fprp(iel,l)  *  (l.OdO  +  fevn(iel,2))  *  (l.OdO  +  fevn(iel,3)) 

&  +  fprp(iel,2)  *  (l.OdO  +  fevn(iel,l))  *  (l.OdO  +  fevn(iel,3)) 

&  +  fprp(iel,3)  *  (l.OdO  +  fevn(iel,l))  *  (l.OdO  +  fevn(iel,2))  ) 

fk(iel)  =  xrhot(iel)  *  cs(iel)**2  * 

&  (  fprp(iel,2)  *  fprp(iel,3)  *  (l.OdO  +  fevn(iel,l)) 

&  +  fprp(iel,l)  *  fprp(iel,3)  *  (l.OdO  +  fevn(iel,2)) 

&  +  fprp(iel,l)  *  fprp(iel,2)  *  (l.OdO  +  fevn(iel,3))  ) 

fl(iel)  =  xrhot(iel)  *  cs(iel)**3  * 

&  fprp(iel,l)  *  fprp(iel,2)*  fprp(iel,3) 

c 

c -  compute  generalized  newmark  time  constants 

c 

tcl(iel)  =  1 . 0d0/(beta*deltat**2) 
tc2(iel)  =  1 . OdO/(beta*deltat) 
tc3(iel)  =  1 . 0d0/(2 . OdO*beta)  -  l.OdO 
tc4(iel)  =  gamma/ (bet a*deltat) 
tc5(iel)  =  (gamma/beta)  -  l.OdO 
tc6(iel)  =  deltat*( (gamma/2. OdO*beta)  -  l.OdO) 
c 

c -  compute  kappa  constants 

c 

kl(iel)  =  fm(iel)*tcl(iel)  +  f c(iel)*tc4(iel) 

k2(iel)  =  fm(iel)*tc2(iel)  +  f c(iel)*tc5(iel) 

k3(iel)  =  fm(iel)*tc3(iel)  +  f c(iel)*tc6(iel) 

c 

6101  continue 
c 

c -  initialize  time  matrices  at  first  time  step  only 

c 


ill 


if  (curstep  .eq.  1  .and.  abs (curstep-lastepl)  .ne.  0)  then 
c 

do  6301  iel  =  ielstrt , ielend 
do  6311  icomp  =  l,nca 
do  6312  jcomp  =  l,nca 
sum_e(iel, icomp, j comp)  =  O.OdO 
sum2_e (iel, icomp, jcomp)  =  O.OdO 
sum_T(iel, icomp, jcomp)  =  O.OdO 
sum2_T (iel, icomp, jcomp)  =  O.OdO 
sum_uij (iel, icomp, jcomp)  =  O.OdO 
etn(iel, icomp, jcomp)  =  O.OdO 
etnpl (iel, icomp, jcomp)  =  O.OdO 

sum_u( iel, icomp, jcomp)  =  O.OdO 
lastepl  =  curstep 
6312  continue 
6311  continue 
6301  continue 


c 


c 

c 

c- 

c 


c 

c- 

c 

c 

c 


c 


call  strain (uxyz,  numel,  ncomp,  nca,  nsol,  deltat ,  ielstrt, 
ielend,  fprp,  fevn,  cs, 

Ecur,  sum_e,  sum2_e,  sum_ui j ,  etn,  etnpl,  curstep) 
else  if  (abs (curstep-lastepl)  .ne.  0)  then 


—  build  strain  vector  after  the  first  time  step 

call  strain (uxyz,  numel,  ncomp,  nca,  nsol,  deltat,  ielstrt, 
ielend,  fprp,  fevn,  cs, 

Ecur,  sum_e,  sum2_e,  sum_ui j ,  etn,  etnpl,  curstep) 


build  stress  vector 


do  6401  iel  =  ielstrt , ielend 


Tcur (iel ,1,1) 
Tcur(iel, 1,2) 
Tcur (iel ,1,3) 
Tcur(iel,2, 1) 
Tcur(iel,2,2) 
Tcur(iel,2,3) 
Tcur(iel,3, 1) 
Tcur(iel,3,2) 
Tcur(iel,3,3) 


xmut (iel) *Ecur (iel ,1,1) 
xmut (iel) *Ecur (iel ,1,2) 
xmut (iel) *Ecur( iel ,1,3) 
xmut(iel)*Ecur(iel,2, 1) 
xmut(iel)*Ecur(iel,2,2) 
xmut (iel) *Ecur (iel ,2,3) 
xmut(iel)*Ecur(iel,3, 1) 
xmut (iel) *Ecur( iel ,3,2) 
xmut (iel) *Ecur( iel ,3,3) 


Tcur (iel , 1 , 1)  =  Tcur (iel , 1 , 1)  +  xlambdat (iel) 
* (Ecur (iel, 1 , 1)+Ecur(iel,2,2)+Ecur(iel,3,3) ) 
Tcur(iel,2,2)  =  Tcur(iel,2,2)  +  xlambdat (iel) 
* (Ecur (iel, 1 , 1)+Ecur(iel,2,2)+Ecur(iel,3,3) ) 
Tcur (iel, 3, 3)  =  Tcur (iel, 3, 3)  +  xlambdat (iel) 
* (Ecur (iel, 1 , 1)+Ecur(iel,2,2)+Ecur(iel,3,3) ) 


c 

c  6401  continue 


c 

c 

c 

c 

c 


sum_T  =  sum_T  +  Tcur 
sum2_T  =  (curstep+curstep**2) *Tcur 
sum_u  =  sum_u  +  u 
lastepl  =  curstep 
endif 


c 

if(idiag  .eq.  GETPRECDIAG  .or.  passnum  .eq.  GETLHS)  then 
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do  8101  iel  =  ielstrt , ielend 


load  up  the  lhs  coefficients 


note  the  coefficients  are  stored  as  follows 
lambda  =  xlambdat (iel) 

mu  =  xmut(iel) 

rho  =  xrhot(iel) 


if  (pml_on  .eq.  1)  then 


a00(iel,l,l)  =  kl(iel)+fk(iel)+0.5*deltat*fl(iel) 
a00(iel,2,2)  =  kl(iel)+fk(iel)+0.5*deltat*fl(iel) 
a00(iel,3,3)  =  kl(iel)+fk(iel)+0.5*deltat*fl(iel) 

all (iel, 1,1)  =  (2.0d0  *  xmut(iel)  +  xlambdat (iel) ) 

*  ((l.OdO  +  fevn(iel,2))*(l .OdO  +  fevn(iel,3)) 
+  0. 5*deltat* (cs (iel) * (fprp(iel ,2) 

*  (l.OdO  +  fevn(iel,3))  +  fprp(iel,3) 

*  (l.OdO  +  fevn(iel,2))))) 

*  ((l.OdO  +  fevn(iel, 1))*(1 .OdO  +  fevn(iel,l)) 

+  0.5*deltat*(cs(iel)*(fprp(iel, 1) 

*  (l.OdO  +  fevn(iel,l))  +  fprp(iel,l) 

*  (l.OdO  +  fevn(iel,l)))))**(-l) 

*  (l.OdO  +  fevn(iel,l)) 

+  0.5*deltat*cs(iel)*fprp(iel,l) 

a22(iel,l,l)  =  xmut(iel) 

*  ((l.OdO  +  fevn(iel,3))*(l .OdO  +  fevn(iel,l)) 
+  0. 5*deltat* (cs (iel) * (fprp(iel ,3) 

*  (l.OdO  +  fevn(iel,l))  +  fprp(iel,l) 

*  (l.OdO  +  f evn(iel ,3) ) ) ) ) 

*  ((l.OdO  +  fevn(iel, 1))*(1 .OdO  +  fevn(iel,2)) 

+  0 . 5*deltat* (cs (iel) * (fprp(iel , 1) 

*  (l.OdO  +  fevn(iel,2))  +  fprp(iel,2) 

*  (l.OdO  +  fevn(iel,l)))))**(-l) 

*  (l.OdO  +  fevn(iel,l)) 

+  0.5*deltat*cs(iel)*fprp(iel,l) 

al2(iel,l,2)  =  xlambdat (iel) 

*  ((l.OdO  +  fevn(iel,2))*(l .OdO  +  fevn(iel,3)) 
+  0 . 5*deltat* (cs (iel) * (fprp(iel ,2) 

*  (l.OdO  +  fevn(iel,3))  +  fprp(iel,3) 

*  (l.OdO  +  fevn(iel,2))))) 

*  ((l.OdO  +  fevn(iel, 1))*(1 .OdO  +  fevn(iel,l)) 

+  0 . 5*deltat* (cs (iel) * (fprp(iel , 1) 

*  (l.OdO  +  fevn(iel,l))  +  fprp(iel,l) 

*  (l.OdO  +  fevn(iel, !)))))**(-!) 


(l.OdO  +  fevn(iel,l)) 

+  0.5*deltat*cs(iel)*fprp(iel,l) 
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a21(iel,l,2)  =  xmut(iel) 


a33(iel ,1,1)  = 


al3(iel ,1,3)  = 


a31 (iel ,1,3)  = 


*  ((l.OdO  +  fevn(iel,3))*(l .OdO  +  fevn(iel,l)) 
+  0 . 5*deltat* (cs (iel) * (fprp(iel ,3) 

*  (l.OdO  +  fevn(iel,l))  +  fprp(iel,l) 

*  (l.OdO  +  f evn(iel ,3) ) ) ) ) 

*  ((l.OdO  +  fevn(iel, 1))*(1 .OdO  +  fevn(iel,2)) 

+  0. 5*deltat* (cs (iel) * (fprp(iel , 1) 

*  (l.OdO  +  fevn(iel,2))  +  fprp(iel,2) 

*  (l.OdO  +  fevn(iel,l)))))**(-l) 

*  (l.OdO  +  fevn(iel,2)) 

+  0.5*deltat*cs(iel)*fprp(iel,2) 

xmut (iel) 

*  ((l.OdO  +  fevn(iel, 1))*(1 .OdO  +  fevn(iel,2)) 
+  0 . 5*deltat* (cs (iel) * (fprp(iel , 1) 

*  (l.OdO  +  fevn(iel,2))  +  fprp(iel,2) 

*  (l.OdO  +  f evn(iel , 1) ) ) ) ) 

*  ((l.OdO  +  fevn(iel, 1))*(1 .OdO  +  fevn(iel,3)) 

+  0. 5*deltat* (cs (iel) * (fprp(iel , 1) 

*  (l.OdO  +  fevn(iel,3))  +  fprp(iel,3) 

*  (l.OdO  +  fevn(iel,l)))))**(-l) 

*  (l.OdO  +  fevn(iel,l)) 

+  0.5*deltat*cs(iel)*fprp(iel,l) 

xlambdat (iel) 

*  ((l.OdO  +  fevn(iel,2))*(l .OdO  +  fevn(iel,3)) 
+  0. 5*deltat* (cs (iel) * (fprp(iel ,2) 

*  (l.OdO  +  fevn(iel,3))  +  fprp(iel,3) 

*  (l.OdO  +  fevn(iel,2))))) 

*  ((l.OdO  +  fevn(iel, 1))*(1 .OdO  +  f evn(iel , 1) ) 

+  0 . 5*deltat* (cs (iel) * (fprp(iel , 1) 

*  (l.OdO  +  fevn(iel,l))  +  fprp(iel,l) 

*  (l.OdO  +  fevn(iel,l)))))**(-l) 

*  (l.OdO  +  fevn(iel,l)) 

+  0.5*deltat*cs(iel)*fprp(iel,l) 

xmut (iel) 

*  ((l.OdO  +  fevn(iel,l))*(1.0d0  +  fevn(iel,2)) 
+  0.5*deltat*(cs(iel)*(fprp(iel, 1) 

*  (l.OdO  +  fevn(iel,2))  +  fprp(iel,2) 

*  (l.OdO  +  f evn(iel , 1) ) ) ) ) 

*  ((l.OdO  +  fevn(iel, 1))*(1 .OdO  +  f evn(iel ,3) ) 

+  0. 5*deltat* (cs (iel) * (fprp(iel , 1) 

*  (l.OdO  +  fevn(iel,3))  +  fprp(iel,3) 

*  (l.OdO  +  fevn(iel,l)))))**(-l) 

*  (l.OdO  +  fevn(iel,3)) 

+  0.5*deltat*cs(iel)*fprp(iel,3) 


all(iel,2,2)  =  xmut(iel) 


a22 (iel ,2,2)  = 


al2(iel,2,l)  = 


a21(iel,2,l)  = 


a33(iel,2,2)  = 


*  ((l.OdO  +  fevn(iel,2))*(l .OdO  +  fevn(iel,3)) 
+  0. 5*deltat* (cs (iel) * (fprp(iel , 2) 

*  (l.OdO  +  fevn(iel,3))  +  fprp(iel,3) 

*  (l.OdO  +  fevn(iel,2))))) 

*  ((l.OdO  +  fevn(iel,2))*(l .OdO  +  fevn(iel,l)) 

+  0. 5*deltat* (cs (iel) * (fprp(iel ,2) 

*  (l.OdO  +  fevn(iel,l))  +  fprp(iel,l) 

*  (l.OdO  +  fevn(iel,2)))))**(-l) 

*  (l.OdO  +  fevn(iel,2)) 

+  0.5*deltat*cs(iel)*fprp(iel,2) 

(2. OdO  *  xmut(iel)  +  xlambdat (iel) ) 

*  ((l.OdO  +  fevn(iel,3))*(l .OdO  +  fevn(iel, 1)) 
+  0. 5*deltat* (cs (iel) * (fprp(iel ,3) 

*  (l.OdO  +  fevn(iel,l))  +  fprp(iel,l) 

*  (l.OdO  +  f evn(iel ,3) ) ) ) ) 

*  ((l.OdO  +  fevn(iel,2))*(l .OdO  +  fevn(iel,2)) 

+  0 . 5*deltat* (cs (iel) * (fprp(iel , 2) 

*  (l.OdO  +  fevn(iel,2))  +  fprp(iel,2) 

*  (l.OdO  +  fevn(iel,2)))))**(-l) 

*  (l.OdO  +  fevn(iel,2)) 

+  0.5*deltat*cs(iel)*fprp(iel,2) 

xmut (iel) 

*  ((l.OdO  +  fevn(iel,2))*(l .OdO  +  fevn(iel,3)) 
+  0. 5*deltat* (cs (iel) * (fprp(iel , 2) 

*  (l.OdO  +  fevn(iel,3))  +  fprp(iel,3) 

*  (l.OdO  +  fevn(iel,2))))) 

*  ((l.OdO  +  fevn(iel,2))*(l .OdO  +  f evn(iel, 1)) 

+  0. 5*deltat* (cs (iel) * (fprp(iel ,2) 

*  (l.OdO  +  fevn(iel,l))  +  fprp(iel,l) 

*  (l.OdO  +  fevn(iel,2)))))**(-l) 

*  (l.OdO  +  fevn(iel,l)) 

+  0.5*deltat*cs(iel)*fprp(iel,l) 

xlambdat (iel) 

*  ((l.OdO  +  fevn(iel,3))*(l .OdO  +  f evn(iel, 1)) 
+  0 . 5*deltat* (cs (iel) * (fprp(iel ,3) 

*  (l.OdO  +  fevn(iel,l))  +  fprp(iel,l) 

*  (l.OdO  +  f evn(iel ,3) ) ) ) ) 

*  ((l.OdO  +  fevn(iel,2))*(l .OdO  +  fevn(iel,2)) 

+  0. 5*deltat* (cs (iel) * (fprp(iel ,2) 

*  (l.OdO  +  fevn(iel,2))  +  fprp(iel,2) 

*  (l.OdO  +  fevn(iel,2)))))**(-l) 

*  (l.OdO  +  fevn(iel,2)) 

+  0. 5*deltat*cs(iel)*fprp(iel ,2) 

xmut (iel) 

*  ((l.OdO  +  fevn(iel, 1))*(1 .OdO  +  fevn(iel,2)) 
+  0 . 5*deltat* (cs (iel) * (fprp(iel , 1) 

*  (l.OdO  +  fevn(iel,2))  +  fprp(iel,2) 

*  (l.OdO  +  f evn(iel , 1) ) ) ) ) 


a23(iel,2,3)  = 


a32(iel,2,3)  = 


al3(iel,3, 1)  = 


a31(iel,3,l)  = 


*  ((l.OdO  +  fevn(iel,2))*(l .OdO  +  fevn(iel,3)) 

+  0. 5*deltat* (cs (iel) * (fprp(iel , 2) 

*  (l.OdO  +  fevn(iel,3))  +  fprp(iel,3) 

*  (l.OdO  +  fevn(iel,2)))))**(-l) 

*  (l.OdO  +  fevn(iel,2)) 

+  0.5*deltat*cs(iel)*fprp(iel,2) 

xlambdat (iel) 

*  ((l.OdO  +  fevn(iel,3))*(l .OdO  +  fevn(iel,l)) 
+  0. 5*deltat* (cs (iel) * (fprp(iel ,3) 

*  (l.OdO  +  fevn(iel,l))  +  fprp(iel,l) 

*  (l.OdO  +  f evn(iel ,3) ) ) ) ) 

*  ((l.OdO  +  fevn(iel,2))*(l .OdO  +  fevn(iel,2)) 

+  0. 5*deltat* (cs (iel) * (fprp(iel , 2) 

*  (l.OdO  +  fevn(iel,2))  +  fprp(iel,2) 

*  (l.OdO  +  fevn(iel,2)))))**(-l) 

*  (l.OdO  +  fevn(iel,2)) 

+  0.5*deltat*cs(iel)*fprp(iel,2) 

xmut (iel) 

*  ((l.OdO  +  fevn(iel, 1))*(1 .OdO  +  fevn(iel,2)) 
+  0.5*deltat*(cs(iel)*(fprp(iel, 1) 

*  (l.OdO  +  fevn(iel,2))  +  fprp(iel,2) 

*  (l.OdO  +  f evn(iel , 1) ) ) ) ) 

*  ((l.OdO  +  fevn(iel,2))*(l .OdO  +  fevn(iel,3)) 

+  0 . 5*deltat* (cs (iel) * (fprp(iel ,2) 

*  (l.OdO  +  fevn(iel,3))  +  fprp(iel,3) 

*  (l.OdO  +  fevn(iel,2)))))**(-l) 

*  (l.OdO  +  fevn(iel,3)) 

+  0.5*deltat*cs(iel)*fprp(iel,3) 


xmut (iel) 

*  ((l.OdO  +  fevn(iel,2))*(l .OdO  +  fevn(iel,3)) 
+  0. 5*deltat* (cs (iel) * (fprp(iel , 2) 

*  (l.OdO  +  fevn(iel,3))  +  fprp(iel,3) 

*  (l.OdO  +  fevn(iel,2))))) 

*  ((l.OdO  +  fevn(iel,3))*(l .OdO  +  fevn(iel,l)) 

+  0. 5*deltat* (cs (iel) * (fprp(iel ,3) 

*  (l.OdO  +  fevn(iel,l))  +  fprp(iel,l) 

*  (l.OdO  +  fevn(iel,3)))))**(-l) 

*  (l.OdO  +  fevn(iel,l)) 

+  0.5*deltat*cs(iel)*fprp(iel,l) 

xlambdat (iel) 

*  ((l.OdO  +  fevn(iel, 1))*(1 .OdO  +  fevn(iel,2)) 
+  0.5*deltat*(cs(iel)*(fprp(iel, 1) 

*  (l.OdO  +  fevn(iel,2))  +  fprp(iel,2) 

*  (l.OdO  +  f evn(iel , 1) ) ) ) ) 

*  ((l.OdO  +  fevn(iel,3))*(l .OdO  +  fevn(iel,3)) 


+  0. 5*deltat* (cs (iel) * (fprp(iel ,3) 


*  (l.OdO  +  f evn(iel ,3) )  +  fprp(iel,3) 

*  (l.OdO  +  fevn(iel,3)))))**(-l) 


a23(iel,3,2)  = 


a32(iel,3,2)  = 


all(iel,3,3)  = 


a22(iel,3,3)  = 


*  (l.OdO  +  fevn(iel,3)) 

+  0. 5*deltat*cs(iel)*fprp(iel ,3) 

xmut (iel) 

*  ((l.OdO  +  fevn(iel,3))*(l .OdO  +  fevn(iel,l)) 
+  0. 5*deltat* (cs (iel) * (fprp(iel ,3) 

*  (l.OdO  +  fevn(iel,l))  +  fprp(iel,l) 

*  (l.OdO  +  f evn(iel ,3) ) ) ) ) 

*  ((l.OdO  +  fevn(iel,3))*(l .OdO  +  f evn(iel , 2) ) 

+  0. 5*deltat* (cs (iel) * (fprp(iel ,3) 

*  (l.OdO  +  fevn(iel,2))  +  fprp(iel,2) 

*  (l.OdO  +  fevn(iel,3)))))**(-l) 

*  (l.OdO  +  fevn(iel,2)) 

+  0.5*deltat*cs(iel)*fprp(iel,2) 

xlambdat (iel) 

*  ((l.OdO  +  fevn(iel, 1))*(1 .OdO  +  f evn(iel , 2) ) 

+  0.5*deltat*(cs(iel)*(fprp(iel, D 

*  (l.OdO  +  fevn(iel,2))  +  fprp(iel,2) 

*  (l.OdO  +  f evn(iel , 1) ) ) ) ) 

*  ((l.OdO  +  fevn(iel,3))*(l .OdO  +  fevn(iel,3)) 

+  0. 5*deltat* (cs (iel) * (fprp(iel ,3) 

*  (l.OdO  +  fevn(iel,3))  +  fprp(iel,3) 

*  (l.OdO  +  fevn(iel,3)))))**(-l) 

*  (l.OdO  +  fevn(iel,3)) 

+  0.5*deltat*cs(iel)*fprp(iel,3) 

xmut (iel) 

*  ((l.OdO  +  fevn(iel,2))*(l .OdO  +  fevn(iel,3)) 
+  0 . 5*deltat* (cs (iel) * (fprp(iel ,2) 

*  (l.OdO  +  fevn(iel,3))  +  fprp(iel,3) 

*  (l.OdO  +  fevn(iel,2))))) 

*  ((l.OdO  +  fevn(iel,3))*(l .OdO  +  fevn(iel,l)) 

+  0. 5*deltat* (cs (iel) * (fprp(iel ,3) 

*  (l.OdO  +  fevn(iel,l))  +  fprp(iel,l) 

*  (l.OdO  +  fevn(iel,3)))))**(-l) 

*  (l.OdO  +  fevn(iel,3)) 

+  0. 5*deltat*cs(iel)*fprp(iel ,3) 

xmut (iel) 

*  ((l.OdO  +  fevn(iel,3))*(l .OdO  +  fevn(iel,l)) 
+  0 . 5*deltat* (cs (iel) * (fprp(iel ,3) 

*  (l.OdO  +  fevn(iel,l))  +  fprp(iel,l) 

*  (l.OdO  +  f evn(iel ,3) ) ) ) ) 

*  ((l.OdO  +  fevn(iel,3))*(l .OdO  +  fevn(iel,2)) 

+  0. 5*deltat* (cs (iel) * (fprp(iel ,3) 

*  (l.OdO  +  fevn(iel,2))  +  fprp(iel,2) 

*  (l.OdO  +  fevn(iel,3)))))**(-l) 


c 

c 

c 


+  0.5*deltat*cs(iel)*fprp(iel,3) 
a33(iel,3,3)  =  (2.0d0  *  xmut(iel)  +  xlambdat (iel) ) 

*  ((l.OdO  +  fevn(iel, 1))*(1 .OdO  +  fevn(iel,2)) 
+  0.5*deltat*(cs(iel)*(fprp(iel, 1) 

*  (l.OdO  +  fevn(iel,2))  +  fprp(iel,2) 

*  (l.OdO  +  f evn(iel , 1) ) ) ) ) 

*  ((l.OdO  +  fevn(iel,3))*(l .OdO  +  fevn(iel,3)) 

+  0. 5*deltat* (cs (iel) * (fprp(iel ,3) 

*  (l.OdO  +  f evn(iel ,3) )  +  fprp(iel,3) 

*  (l.OdO  +  fevn(iel,3)))))**(-l) 

*  (l.OdO  +  fevn(iel,3)) 

+  0.5*deltat*cs(iel)*fprp(iel,3) 


else 

c 

aOO(iel,l,l)  =  xrhot (iel) *tcl (iel) 
a00(iel,2,2)  =  xrhot (iel) *tcl (iel) 
a00(iel,3,3)  =  xrhot (iel) *tcl (iel) 
c 

all (iel, 1,1)  =  2. OdO  *  xmut(iel)  +  xlambdat (iel) 

a22(iel,l,l)  =  xmut(iel) 
al2(iel,l,2)  =  xlambdat (iel) 
a21(iel,l,2)  =  xmut(iel) 
a33(iel,l,l)  =  xmut(iel) 
al3(iel,l,3)  =  xlambdat (iel) 
a31(iel,l,3)  =  xmut(iel) 
c 

all(iel,2,2)  =  xmut(iel) 

a22(iel,2,2)  =  2. OdO  *  xmut(iel)  +  xlambdat (iel) 

al2(iel,2,l)  =  xmut(iel) 
a21(iel,2,l)  =  xlambdat (iel) 
a33(iel,2,2)  =  xmut(iel) 
a23(iel,2,3)  =  xlambdat (iel) 
a32(iel,2,3)  =  xmut(iel) 
c 

al3(iel,3,l)  =  xmut(iel) 
a31(iel,3,l)  =  xlambdat (iel) 
a23(iel,3,2)  =  xmut(iel) 
a32(iel,3,2)  =  xlambdat (iel) 
all(iel,3,3)  =  xmut(iel) 
a22(iel,3,3)  =  xmut(iel) 

a33(iel,3,3)  =  2. OdO  *  xmut(iel)  +  xlambdat (iel) 

c 

endif 

c 

8101  continue 
c 

endif 

c 

c 

c -  load  up  the  rhs  coefficients 

c 

if(iresid  .eq.  GETRHSORRESIDUAL  .or.  passnum  .eq.  GETLHS)  then 
c 

if  (pml_on  .eq.  1)  then 
c 

do  9101  iel  =  ielstrt , ielend 
do  9111  icomp  =  l,nca 
c 

f(iel, icomp)  =  (kl(iel)  *  u(iel, icomp, soltn) 
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+  k2(iel)  *  u(iel,icomp,veltn) 

+  k3(iel)  *  u(iel,icomp,acctn) 

-  f 1 (iel) *deltat*sum_u(iel , icomp , soltn) ) 

continue 

fx(iel,l)  =  -(((l.OdO  +  fevn(iel,2))*(l .OdO  +  fevn(iel,3)) 
+  0 . 5*deltat* (cs (iel) * (f prp (iel , 2) 

*  (l.OdO  +  fevn(iel,3))  +  fprp(iel,3) 

*  (l.OdO  +  fevn(iel,2))))) 

*Tcur(iel, 1 , 1) 

+  (l.OdO  +  fevn(iel,2))*(1.0d0  +  fevn(iel,3)) 
*deltat*sum_T(iel ,1,1) 

+  ( (cs (iel)**2) *fprp(iel , 2) *fprp(iel ,3) ) 

* (deltat**2) *sum2_T(iel ,1,1)) 

fx(iel,2)  =  -(((l.OdO  +  fevn(iel,3))*(l .OdO  +  fevn(iel,l)) 
+  0 . 5*deltat* (cs (iel) * (fprp(iel ,3) 

*  (l.OdO  +  fevn(iel,l))  +  fprp(iel,l) 

*  (l.OdO  +  f evn(iel ,3) ) ) ) ) 

*Tcur(iel, 1 ,2) 

+  (l.OdO  +  f evn(iel,3) ) *(1 .OdO  +  fevn(iel,l)) 
*deltat*sum_T(iel ,1,2) 

+  ((cs(iel)**3)*fprp(iel,2)*fprp(iel, 1)) 

* (deltat**2) *sum2_T(iel ,1,2)) 

fx(iel,3)  =  -(((l.OdO  +  fevn(iel,l))*(1.0d0  +  fevn(iel,2)) 
+  0 . 5*deltat* (cs (iel) * (fprp(iel , 1) 

*  (l.OdO  +  fevn(iel,2))  +  fprp(iel,2) 

*  (l.OdO  +  f evn(iel, 1) ))) ) 

*Tcur (iel ,1,3) 

+  (l.OdO  +  fevn(iel,l))*(1.0d0  +  fevn(iel,2)) 
*deltat*sum_T(iel, 1,3) 

+  ( (cs (iel)**l) *fprp(iel ,2) *fprp(iel ,2) ) 

* (deltat**2) *sum2_T(iel , 1,3)) 

print  * , ’fx* ,fx 

fy(iel,l)  =  -(((l.OdO  +  fevn(iel,2))*(l .OdO  +  fevn(iel,3)) 
+  0 . 5*deltat* (cs (iel) * (fprp(iel ,2) 

*  (l.OdO  +  fevn(iel,3))  +  f prp (iel, 3) 

*  (l.OdO  +  fevn(iel,2))))) 

*Tcur(iel,2,l) 

+  (l.OdO  +  fevn(iel,2))*(l .OdO  +  fevn(iel,3)) 
*deltat*sum_T(iel ,2,1) 

+  ( (cs (iel)**2) *fprp(iel ,2) *fprp(iel ,3) ) 

* (deltat**2) *sum2_T(iel ,2,1)) 

fy(iel,2)  =  -(((l.OdO  +  fevn(iel,3))*(l .OdO  +  fevn(iel,l)) 
+  0 . 5*deltat* (cs (iel) * (fprp(iel ,3) 

*  (l.OdO  +  fevn(iel,l))  +  fprp(iel,l) 

*  (l.OdO  +  f evn(iel ,3) ) ) ) ) 

*Tcur(iel,2,2) 

+  (l.OdO  +  fevn(iel,3))*(l .OdO  +  fevn(iel,l)) 
*deltat*sum_T(iel ,2,2) 

+  ((cs(iel)**3)*fprp(iel,2)*fprp(iel, 1)) 

* (deltat**2) *sum2_T(iel ,2,2)) 

fy(iel,3)  =  -(((l.OdO  +  fevn(iel,l))*(1.0d0  +  fevn(iel,2)) 
+  0 . 5*deltat* (cs (iel) * (fprp(iel , 1) 

*  (l.OdO  +  fevn(iel,2))  +  fprp(iel,2) 

*  (l.OdO  +  f evn(iel , 1) ) ) ) ) 

*Tcur(iel,2,3) 

+  (l.OdO  +  fevn(iel, 1))*(1 .OdO  +  fevn(iel,2)) 
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c 


c 


c 

c 

c 

9101 


*deltat*sum_T (iel ,2,3) 

+  ( (cs(iel)**l) *fprp(iel,2) *fprp(iel,2) ) 
* (deltat**2) *sum2_T(iel ,2,3)) 


print  * , ’  fy’ ,fy 

fz(iel,l)  =  -(((l.OdO  +  fevn(iel,2))*(l .OdO  +  fevn(iel,3)) 
+  0 . 5*deltat* (cs(iel) * (fprp(iel ,2) 

*  (l.OdO  +  fevn(iel,3))  +  fprp(iel,3) 

*  (l.OdO  +  f evn(iel , 2) ) ) ) ) 

*Tcur(iel,3, 1) 

+  (l.OdO  +  fevn(iel,2))*(1.0d0  +  fevn(iel,3)) 
*deltat*sum_T(iel,3, 1) 

+  ( (cs (iel) **2) *fprp(iel ,2) *fprp(iel ,3) ) 

* (deltat**2) *sum2_T(iel ,3,1)) 

fz(iel,2)  =  -(((l.OdO  +  fevn(iel,3))*(l .OdO  +  fevn(iel,l)) 
+  0 . 5*deltat* (cs (iel) * (fprp(iel ,3) 

*  (l.OdO  +  fevn(iel,l))  +  fprp(iel,l) 

*  (l.OdO  +  f evn(iel ,3) ) ) ) ) 

*Tcur(iel,3,2) 

+  (l.OdO  +  fevn(iel,3))*(l .OdO  +  fevn(iel,l)) 
*deltat*sum_T (iel ,3,2) 

+  ( (cs (iel) **3) *fprp(iel , 2) *fprp(iel , 1) ) 

* (deltat**2) *sum2_T(iel ,3,2)) 


fz(iel,3)  =  -(((l.OdO  +  fevn(iel,l))*(1.0d0  +  fevn(iel,2)) 
+  0 . 5*deltat* (cs (iel) * (fprp(iel , 1) 

*  (l.OdO  +  fevn(iel,2))  +  fprp(iel,2) 

*  (l.OdO  +  f evn(iel, 1) ))) ) 

*Tcur(iel,3,3) 

+  (l.OdO  +  fevn(iel,l))*(1.0d0  +  fevn(iel,2)) 
*deltat*sum_T (iel ,3,3) 

+  ( (cs (iel) **1) *fprp(iel ,2) *fprp(iel ,2) ) 

* (deltat**2) *sum2_T(iel ,3,3)) 


print  *, ’fz} ,fz 
continue 


c 

else 

c 

do  9151  icomp  =  l,nca 

do  9161  iel  =  ielstrt , ielend 
c 

f (iel, icomp)  =  xrhot(iel)  *  ( 

&  u(iel , icomp, soltn)  /  (beta  *  deltat**2)  + 

&  u(iel , icomp, veltn)  /  (beta  *  deltat)  + 

&  u(iel , icomp, acctn)  *  (1 . 0/(2 . 0*beta)  -  1.0)  ) 

c 

fx(iel , icomp)  =  O.OdO 
fy (iel , icomp)  =  O.OdO 
fz (iel , icomp)  =  O.OdO 
c 

9161  continue 
9151  continue 
c 

endif 

endif 

c 

c 

return 

end 

c 
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strain 


c 

subroutine  strain (uxyz,  numel,  ncomp,  nca,  nsol,  deltat,  ielstrt, 
ielend,  fprp,  fevn,  cs, 

Ecur,  sum_e,  sum2_e,  sum_ui j ,  etn,  etnpl,  curstep) 
c 

(^3f:3|e^c3f;3|c^c3f;^e^c^e^e^c3|e^e^o|c^e^e3|e^c^e3|c^c3f;3|e^c^e^e^c3|e3|e^c3|c^e3f;3|e^e^e3|e^c3f;3|e^c^e3|c^c3f;^e^c3|e^e3f;3|e^e^o|c^c3f;3|e^c^e^e^c3f;3|e^c3|c^e^c3|e^e 


C*  * 

c*  Routine:  strain  * 

c*  Purpose:  computes  the  strain  solution,  velocity  and  * 

c*  acceleration  vectors  at  time  T(n+1)  based  on  the  * 

c*  solution,  velocity,  and  acceleration  at  time  T(n).  * 

c*  * 

c*  Parameters:  * 

c*  I  uxyz  -  * 

c*  I  Ecur  -  current  strain  * 

c*  I  etn  -  strain  at  time  Tn  * 

c*  I  etnpl  -  strain  at  time  Tn+1  * 

c*  * 

c*  * 


Q^o|e^c3)c^e^c3f;3|e^c3|e^e^c3|e^e3f;3|e^e3f;3|e^c3f;^e3|c3f;3|e^c3|c^e^c3|e^e3f;3|e3|e3f;3|e^c^e3|e^c3f;3|e^c^e3|c^c3|e^e^c3|e^e3f;3|e^c3f;3|c^c3f;3|e^c^e^e^c3|e3|e^c3|c^e3f;3|o|e 

c 

#include  "PH-std.blk" 
c 
c 
c 

INTTYPE  numel,  ncomp,  nca,  nsol,  iel,  ielstrt,  ielend 
INTTYPE  icomp,  jcomp,  k,  soltn,  veltn,  acctn,  curstep 
c 

REALTYPE  deltat 

REALTYPE  uxyz (numel, ncomp, nsol, 3) , 

&  Ecur (MAXELEMBATCH , ncomp , ncomp) , 

&  etn (MAXELEMBATCH, ncomp, ncomp) , 

&  etnpl (MAXELEMBATCH, ncomp, ncomp) , 

&  sum_e (MAXELEMBATCH , ncomp , ncomp) , 

&  sum2_e (MAXELEMBATCH , ncomp , ncomp) , 

&  sum_ui j (MAXELEMBATCH , ncomp , ncomp) , 

&  fprp (MAXELEMBATCH , ncomp) , 

&  f evn (MAXELEMBATCH, ncomp) , 

&  cs (MAXELEMBATCH) 

c 

c  print  *, ’inside  strain  subroutine. . . . ’ 

c 

c - 

c 

call  PH_UG_GET_SLOT_ (S0LUTI0N_TN ,  soltn) 
call  PH_UG_GET_SLOT_ (VELOCITY_TN ,  veltn) 
call  PH_UG_GET_SLOT_ (ACCELERATION_TN ,  acctn) 
c 
c 

etn  =  Ecur 
sum_e  =  sum_e  +  etn 
c 

sum2_e  =  (curstep  +  curstep**2) *etn 
c 

do  3001  iel  =  ielstrt , ielend 


c 

c -  building  sum  vector 

c 

sum_uij (iel , 1 , 1)  =  (1 . 0d0*sum_uij (iel , 1 , 1)  + 
&  (  fprp(iel , 1) *uxyz(iel , 1 , soltn, 1) 

&  +  fprp(iel , 1) *uxyz(iel , 1 , soltn, 1) )  ) 

c 

sum_uij (iel , 1 , 2)  =  (1 . 0d0*sum_uij (iel , 1 ,2)  + 
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&  (  fprp(iel , 1) *uxyz(iel , 1 , soltn,2) 

&  +  fprp(iel,2)*uxyz(iel,2,soltn,l))  ) 

c 

sum_uij (iel , 1 ,3)  =  (1 . OdO*sum_uij (iel , 1 ,3)  + 

&  (  fprp(iel , 1) *uxyz(iel , 1 , soltn,3) 

&  +  fprp(iel,3)*uxyz(iel,3,soltn,l))  ) 

c 

sum_uij (iel ,2, 1)  =  (1 . OdO*sum_ui j (iel , 2 , 1)  + 

&  (  fprp(iel,2)*uxyz(iel,2,soltn, 1) 

&  +  fprp(iel , 1) *uxyz(iel , 1 , soltn,2) )  ) 

c 

sum_uij (iel, 2, 2)  =  (1 . OdO*sum_uij (iel, 2, 2)  + 

&  (  fprp(iel,2)*uxyz(iel,2,soltn, 1) 

&  +  fprp(iel,2)*uxyz(iel,2,soltn,2))  ) 

c 

sum_uij (iel ,2,3)  =  (1 . OdO*sum_ui j (iel, 2, 3)  + 

&  (  fprp(iel,2)*uxyz(iel,2,soltn,3) 

&  +  fprp(iel,3)*uxyz(iel,3,soltn,2))  ) 

c 

sum_uij (iel ,3, 1)  =  (1 . OdO*sum_ui j (iel ,3 , 1)  + 

&  (  fprp(iel,3)*uxyz(iel,3,soltn, 1) 

&  +  fprp(iel , 1) *uxyz(iel , 1 , soltn,3) )  ) 

c 

sum_uij (iel ,3,2)  =  (1 . OdO*sum_uij (iel, 3, 2)  + 

&  (  fprp(iel,3)*uxyz(iel,3,soltn,2) 

&  +  fprp(iel,2)*uxyz(iel,2,soltn,3))  ) 

c 

sum_uij (iel ,3,3)  =  (1 . OdO*sum_uij (iel, 3, 3)  + 

&  (  fprp(iel,3)*uxyz(iel,3,soltn,3) 

&  +  fprp(iel,3)*uxyz(iel,3,soltn,3))  ) 

c 

c -  build  strain  vector  for  t(n+l) 

c 

do  9201  icomp  =  l,nca 
do  9211  jcomp  =  l,nca 
c 

etnpl (iel, icomp, jcomp)  =  0.5d0 
&  *((1.0d0  +  f evn(iel, icomp) )*(1 .OdO  +  f evn(iel , j comp) ) 

&  +  0.5*deltat*(cs(iel)*(fprp(iel,icomp)*(1.0d0  +  fevn(iel, jcomp)) 

&  +  fprp(iel, jcomp)*(1.0d0  +  f evn(iel , icomp) ))))** (-1) 

&  *  (1 .0d0*cs (iel) *deltat*sum_uij (iel, icomp, jcomp) 

&  -  2 . 0d0*deltat* (cs (iel) * (fprp(iel , icomp) * (1 . OdO  +  fevn(iel, jcomp) ) 

&  +  fprp(iel , jcomp) * (1 . OdO  +  f evn(iel , icomp) )) ) 

&  *  sum_e (iel , icomp, jcomp) 

&  -  2 . 0d0*deltat**2* ( (cs (iel) **2) *fprp(iel , icomp) *fprp(iel , jcomp) ) 

&  *  sum2_e (iel, icomp, jcomp)  ) 

c 

9211  continue 
9201  continue 
c 

3001  continue 
c 

Ecur  =  etnpl 
c  print  * ,  Ecur 
c 
c 

return 

end 
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APPENDIX.  F  (PROPHLEX  C++  MODULE) 


/*  -  application. c  -  Fri  Jul  26  14:29:13  CDT  1996  - 

*  Copyright:  Computational  Mechanics,  Co.,  Inc.  1992-1996 

*  MAJ  Anthony  N.  Johnson,  Naval  Postgraduate  School 

*  Monterey,  California  93940 


3|e^e3|c^c3|c3|e^c3f;3|e^:3f;^e^c3|e^e^e3|e3|e^e3|e^c3|e3|c^c^e3|e^c^e3|c^c3|e^e^c3|e3|e^e3|e^e3f;3|e^c^e3|e^c4;3|e^c3f;3|e^c3|e^e^e3|e^e^e3|e^c^e3|e^c^ 

/*  #####################################################################  */ 
/*  File:  application. c 
*  ========================  */ 

#include  "PH-std.h" 

#include  <string.h> 


/* 

/*  Function  Prototypes 


*/ 


*/ 


#include  "PH-util.h" 
#include  "PH-register .h" 
#include  "PH-parameters .h" 
#include  "PH-adapt.h" 


#def ine  APPLICATION 
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/*  the  application  number  */ 


#include  "PH-Fmacros .h" 

#include  "PH-post.h" 

#include  "PH-gui.h" 

#include  "PH-tcl.h" 

#include  "PH-io.h" 

#include  "PH-results-save .h" 

#include  "_release_name .h" 

#include  "application-parameters .h" 
#include  "application-macros .h" 
#include  "application-main .h" 

#include  "application-post .h" 

#include  "application-FtoC .h" 

#include  "applicationFile-Register .h" 


/* 

/*  Global  variable 


*/ 


*/ 


extern  Tcl_Interp  *  tel; 


/*  Internal  Function  Prototypes 

* - 3|e  / 

static  void  application_tcl (void) ; 

static  int  cmdApplication(ClientData  clientData,  Tcl_Interp  *  interp, 
int  argc ,  char  *argv [] ) ; 


/*  Static  Variables 
* - */ 

/*  define  postprocessor  components  and 
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*  type  of  data  needed  for  evaluation 
*/ 


static  struct  PostComp  postApplication []  = 
{ 


{XCOMP,  "X", 

NEED_S0L0NLY>, 

{YCOMP,  "Y", 

NEED_S0L0NLY>, 

{ZCOMP,  "Z", 

NEED.SOLONLY}, 

{PRIMARY 1 , 

"X  Displacements 

",  NEED.SOLONLY} , 

{PRIMARY2 , 

"Y  Displacements 

",  NEED.SOLONLY}, 

{PRIMARY3 , 

"Z  Displacements 

",  NEED.SOLONLY}, 

{SECOND ARY1 , 

"Secondary  1", 

NEED.DERIY} , 

{SECOND ARY2, 

"Secondary  2", 

NEED.DERIY} , 

{SECOND  ARY3, 

"Secondary  3", 

NEED.DERIY} , 

{NEWCOMP 

"New  Comp"  , 

NEED.DERIV} , 

{NEWCOMPWAVE 

,  "Wave  Comp"  , 

NEED.DERIV} , 

{ESTIMATE,  "] 

Error  Indicator", 

ELEMENT.COMP} , 

/*  the  following  is  always  last  */ 

{0,  >\0\  0} 


/*  ##################################################################  */ 
/*  declaration  of  structures  for  HM  output  file  formats  */ 

static  int  Num_Results_Fmts  =  1  ; 

static  ResultsFileFmt  Sav_Results_Fmts [] = 

{ 

{"res",  "HM",  "*.res",  1, 

FM_BINARY, 

PH_WriteHM_Results , 

PH_WriteHM_Mesh, 

PH_WriteHM.AH, 

PH.AppendHM.Step, 

PH.FinishHM.Steps , 

9, 


{ 


{ 

XCOMP 

RF.SCALAR 

1, 

"X">, 

{ 

YCOMP 

RF.SCALAR 

2, 

"Y">, 

{ 

ZCOMP 

RF.SCALAR 

3, 

"Z">, 

{ 

PRIMARY 1 

RF.VECTOR 

4, 

"Displacements"} 

{ 

PRIMARY2 

RF.VECTOR 

4, 

"Displacements"} 

{ 

PRIMARY3 

RF.VECTOR 

4, 

"Displacements"} 

{ 

SECONDARY 1  , 

RF.SCALAR  , 

5, 

"Secondary  1"}, 

{ 

SECOND ARY2  , 

RF.SCALAR  , 

6, 

"Secondary  2"}, 

{ 

SECOND ARY3  , 

RF.SCALAR  , 

7, 

"Secondary  3"}, 

{ 

NEWCOMP 

RF.SCALAR  , 

8, 

"New  Comp"  }, 

{ 

NEWCOMPWAVE  , 

RF.SCALAR  , 

9, 

"Wave  Comp"  } 

/* 

Vector  components  are  also  possible 

{ 

SOLI 

RF.VECTOR 

4, 

"Displacements"} 

{ 

S0L2 

RF.VECTOR 

4, 

"Displacements"} 

124 


{  S0L3 


RF_VECTOR 


4,  "Displacements"}, 


*/ 

} 

} 

}; 

/*  ##################################################################  */ 


/*  Define  the  solution  algorithm  name  and  potential  linear  solvers 

*  to  call.  Note,  the  alorithm  integer  reference  number  for 

*  ALG_TEMPLATE  is  stored  under  KP_ALGORITHM  in  the  parameter  set 

*  KERNEL_PARAMS .  This  parameter  may  be  used  to  branch  between 

*  algorithms  in  a  given  application  e.g.  (linear,  nonlinear,  transient  ...) 
*/ 

static  Alg_Data_Type  algApplication []  = 

{ 

{ALG_TEMPLATE,  "Your  algorithm", 

{SPARSE_OPT,  SPARSE_OPT,  SPARSE_OPT,  SPARSE_OPT,  SPARSE_OPT}} , 

{0,  >\0>} 

}; 


/*  define  the  sequence  of  buttons  and  names 
*  to  appear  on  the  viewport  window 
*/ 

static  struct  PostTopbar  post_menu[]  =  { 


{C0NTR0L_F0RM , 

"Control"} , 

{C0L0R_EDIT_F0RM , 

"C0L0R2" ,  "C0L0R2"} 

{SET_SELECT_FORM , 

"SETS",  "SETS"}, 

{VIEW_F0RM, 

"EYE2"} , 

{C0MP_SELECT0R_F0RM , 

" LoadRe  suit  s " } , 

{SLICE_F0RM, 

"Slice"}, 

{IS0SURF_F0RM, 

"Iso"}, 

{XY_F0RM, 

"XYplot "} , 

{PR0BE_F0RM, 

"Probe"}, 

{L0ADS_PL0T_F0RM , 

"BCs"} , 

{HARDC0PY_F0RM , 

"Hardcopy"}, 

{DEFRM_F0RM , 

"Deform"} , 

{VELVEC_F0RM, 

"VelVec"} , 

{TRACES_F0RM, 

"Traces"} , 

{RESET_ACTI0N, 

"ResetView"} , 

{0,  NULL}, 

}; 

/*#########################################################################  */ 
/*  EXPORTABLE  FUNCTIONS  */ 

/*  #######################################################################  */ 
int  main  (  int  argc,  char  **  argv  ) 

{ 

return  PH_PhlexMain  (  argc ,  argv  ) ; 

} 

/*  ##################################################################  */ 
void  PH_InitializeApplicationDB (void) 

{ 


/* 

*  This  routine  is  the  main  registration  and  initialization  routine 

*  called  by  default  from  the  kernel 
*/ 
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/*  set  the  icon  name  to  appear  on  the  main  form  and 

*  the  release  name 
*/ 

PH_Gui_SetMainForm  (  "ProPHLEX" ,  NULL  ); 

PH_Gui_SetReleaseName (_RELEASE_NAME) ; 

applicationFile_RegisterObjDB  () ;  /*  attach  user  objects  to  database  */ 

/*  Register  the  following; 

*  application  name,  initialization  routine 

*  postprocessor  components  to  evaluate 

*  HM  interface  options 

*  linear  equation  solver  to  link 

*  tel  initialization  routine  for  special  application  parameters 
*/ 

PH_RegisterProblem (APPLICATION,  "Template" , 

SCOPE,  (void  *)  &algApplication) ; 

PH_RegisterPostcomp (APPLICATION,  (void  *)  &post Application,  0,  PRIMARY1, 

NULL  ,  NULL) ; 

PH_RegisterResultsFmts  (APPLICATION,  Num_Results_Fmts ,  Sav_Results_Fmts) ; 
PH_RegisterSolver (  SPARSE_0PT,  "Sparse"  ,  PH_INIT_SPARSE_  );  /*  sparse  solver  */ 
PH_RegisterTclInitialization( (FUNPTR)  application_tcl) ; 

/*  Register  the  following: 

*  the  first  slot  to  use  for  velocity  vectors  in  post 

*  note  for  deformed  plots  the  first  three  slots  are  used  by  default 
*/ 


PH_Gui_SetPostOption  (P0ST_VECT0R_C0MP_ID ,  PRIMARYl) ; 


f pr intf (PH_dbgout ,  " \n\n" ) ; 

f printf (PH_dbgout ,  ">>>>  ====================================  <<<<\n" ) 

f printf  (PH_dbgout ,  "»»  WELCOME  TO  SAFE-T  ««\n") 

f  printf  (PH_dbgout ,  "»»  hp-ADAPTIVE  FINITE  ELEMENT  ANALYSIS  ««\n") 
f printf (PH_dbgout ,  ">>»  TOOL  using  the  PHLEX  kernel  <<<<\n") 

f printf (PH_dbgout ,  ">>>>  ====================================  <<<<\n") 

f printf (PH_dbgout ,  " \n\n" ) ; 


/*  ####################################################################  */ 
void  PH_InitializeApplicationParams (  void  ) 

{ 

/*  this  routine  initializes  and  set  default  values  for; 

*  the  aplication  parameters 

*  the  material  data  base  parameters 

*  the  entries  to  appear  on  the  veiwport  window 
*/ 

SetDef aultApplicationParams () ; 

SetDef  aultMaterialDBO  ; 

SetDef  aultBoundaryCondit ionsDB ( ) ; 

PH_Gui_SetPostMainMenu  (  post_menu  ) ; 

/*  turn  on  initial  mesh  conditioning  */ 
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PH_Gui_ConditionMesh("on" , (double) 0) ; 


> 

/*  ####################################################################  */ 
static  void  application_tcl (void) 

{ 

/* 

*  initialization  of  application  specific  tel  options 

*  in  particular,  the  function  cmdApplication  below  is  registered 

*  as  a  call  back  function  from  the  kernel 
*/ 


Tcl_CreateCommand(tcl,  "command",  cmdApplication,  (ClientData)  NULL, 
(Tcl_CmdDeleteProc  *)  NULL) ; 

} 

/*  ####################################################################  */ 
static  int  cmdApplication(  ClientData  clientData, 

Tcl_Interp  *interp, 
int  argc , 
char  *argv  []  ) 

{ 

/* 

*  application  specific  tel  options  are  defined  in  this  routine 
*/ 


char  ActionType [MAX_PARAM_NAME_LEN] ; 
int  action; 

if  (argc  !=  3)  { 

Tcl_AppendResult (interp,  "wrong  #  args :  should  be  \"",  argv[0] , 
"  action_type  well_number\" " ,  (char  *)  NULL); 
return  TCL_ERR0R; 

> 

strcpy (ActionType ,  argv[l]); 

if  (strcmp (ActionType,  "token")  ==  STRING_MATCH)  { 
action  =  YES; 

}  else  { 

Tcl_AppendResult (interp,  "unable  to  identify  action  type", 

(char  *)  NULL); 
return  TCL_ERR0R; 

> 

return  TCL_0K; 


} 
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APPENDIX.  G  (DISCRETE  CONVOLUTION 

USING  MATLAB) 


clear  M; clear  a; clear  f ; clear  h; clear  p;  close  all 

a=0 . 005; 

j=0; 

etime=2; 

step=0 . 01 ; 

time=0 : step : etime ; 

length(time) ; 

for  s  =  1 : length (time) 

j  =  j+i; 

7.  p(j)  =  -0.055; 

7o  p(j)  =  -2*sourceGaus (time (s) ) ; 
p(j)  =  -2 . 5*derGaus (time (s) ) ; 

end 

r=l; 

length (p) ; 

M=zeros (length(p) ) ; 
length(p) ; 

[m,n] =size (M) ; 
for  i  =  l:m 

for  j  =  l:n 

°/„  M(i,  j)=Heavi(i-j)-Heavi(i-j-l) ; 

h(j)  =  wtime(time(j) ,r) ; 

M(i, j)=wtime(time(j)-time(i) ,r)-wtime(time(j)-time(i)-step,r) ; 

end 

end 

f=p*M; 

subplot (3,1,1) ;plot (time ,f ) 
grid 

title ( ’ Surf ace  Response  lm  away’) 
ylabel ( ’ Displacement ’ ) 

subplot (3,1,2) ;plot (time ,p, ’r ’ ) ; 
grid 

title ( ’ Surf ace  Source') 
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ylabel ( ’ Displacement ’ ) 


subplot (3,1,3) ;plot (time ,h, ’g’ ) ; 
grid 

title( ’Pekeris  1955’) 
xlabel(’time  (seconds)’) 
ylabel ( ’Displacement ’ ) 


f igure(2) 

plot(time,f) 

grid 

title (’ Surf ace  Response  lm  away’) 
ylabel ( ’Displacement ’ ) 
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